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SECTION  I 
INTRODUCTION 


The  optimum  performance  of  a  turbo -propuls ion  system  is  usually 
achieved  when  the  compressor  is  operating  near  its  maximum  pressure  ratio. 

However,  this  optimum  is  generally  not  attainable  because  it  occurs  close  to 
compressor  stall  and  unstable  flow  conditions.  In  actual  operation,  a  stall 
margin  must  be  provided  to  prevent  the  compressor  from  penetrating  the  stall 
boundary  and  developing  destructive  unsteady  flow  phenomena  such  as  rotating 
stall  and  surge.  This  is  usually  done  by  prescheduling  the  engine  controls. 

When  an  aircraft  has  a  varied  flight  envelope,  the  prescheduling  approach  can 
lead  to  the  requirement  for  a  large  stall  margin  to  keep  the  engine  out  of 
stall  under  all  possible  transient  and  steady  flight  conditions.  This  stall 
margin  represents  a  significant  performance  penalty.  Also,  in  many  instances 
of  engine  failure,  rotating  stall  has  been  identified  as  a  precursor  to  destructive 
unsteady  flows  in  an  engine.  Furthermore,  blade  fatigue  considerations  will 
not  allow  a  compressor  to  operate  for  prolonged  periods  in  a  large  amplitude 
rotating  stall  mode.  Clearly  then,  it  is  desirable  to  develop  methods  of 
estimating  the  stall  boundaries  of  a  compressor  and  if  possible  to  develop  an 
engine  control  system  that  can  sense  incipient  destructive  unsteady  flows  in  a 
compressor  and  take  corrective  action  to  prevent  compressor  stall.  Recognition 
of  these  goals  has  been  the  motivation  for  a  continuing  program  of  research  that 
the  Aeropropulsion  Laboratory  has  sponsored  at  Calspan  since  1962.  The  last 
program  at  Calspan  was  carried  out  under  Contract  No.  F33615-76-C-2092  and  the 
results  are  reported  in  Refs.  1  through  3. 

The  work  at  Calspan  has  been  both  theoretical  and  experimental  in 
nature  and  has  been  aimed  at  obtaining  a  sufficient  understanding  of  the  rotating 
stall  phenomenon  such  that  its  onset  and  its  properties  can  be  predicted  and 
controlled.  Demonstrated  progress,  has  been  made  toward  these  goals  in  that  a 
theory  has  been  developed  which  is  capable  of  predicting  inception  of  rotating 
stall  on  a  high  hub-to-tip  ratio  compressor  stage  (rotor  plus  stator)  in  low 
speed  flows  (provided  that  the  appropriate  steady  state  blade  row  performance 
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data  are  available).  In  addition,  a  prototype  rotating  stall  control  system  has 
been  designed  and  demonstrated  successfully  by  tests  conducted  by  Calspan  on  a 
J85-5  turbojet  engine  and  by  the  Air  Force  Aero  Propulsion  Laboratory  on  a  J85-13 
turbojet  engine. 

The  latest  three-year  program  has  the  objectives  of  (1)  develop  an 
analysis  for  a  three-dimensional  time-variant  rotating  stall  and  separation 
theory,  (2)  develop  analyses  for  post-stall  operation/recovery  and  aerodynamical ly 
induced  exotic  metal  combustion  and  (3)  consider  the  effects  of  distortion,  water 
ingestion,  and  nuclear  blasts  on  axial  flow  compressors.  The  work  done  towards 
accomplishing  objective  (1)  is  reported  on  herein.  The  work  done  towards  the 
accomplishments  of  the  remaining  objectives  are  reported  in  Volume  2  of  this 
report . 


The  theoretical  work  performed  under  the  present  program  was  aimed  at 
developing  methods  of  analysis  to  predict  the  unsteady,  viscous,  three- 
dimensional  flow  through  a  blade  row.  The  effects  of  viscosity  and  unsteadiness 
were  considered  separately.  Thus -the  attempt  was  to  develop  a  three-dimensional, 
steady-flow  theory  to  predict  the  turning  and  loss  performance  of  a  blade  row 
and  a  method  to  estimate  the  lossless  unsteady,  three-dimensional  flow  through 
a  blade  row  due  to  inlet  distortion. 

Past  correlations  made  with  the  small  disturbance  stability  theory  in 
incompressible  flows  have  shown  that  if  the  steady-state  loss  and  turning 
performance  of  a  blade  row  or  stage  are  known,  the  rotating  stall  inception 
boundary  could  be  predicted  quite  accurately.  Moreover,  the  two-dimensional 
stability  theory  worked  well  for  distinctly  three-dimensional  flows  if  radially 
averaged  blade  row  performance  was  used.  Consequently,  the  theoretical  work 
under  this  program  concentrated  on  extending  the  stability  theory  to  compressible 
flow  and  initiating  the  development  of  methods  to  predict  the  loss  and  turning 
performance  of  viscous  flow  through  blade  rows.  A  version  of  the  stability 
theory  applicable  to  a  single  blade  row  in  two-dimensional  compressible  flow 
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had  been  previously  formulated  ,  but  no  results  were  presented.  This  analysis 
was  extended  to  consider  two  blade  rows  and  was  implemented  into  a  computer  code. 
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The  basic  philosophy  that  was  chosen  to  analyze  the  viscous  flow  through 
a  blade  row  was  to  split  the  flow  into  an  inviscid  region  and  a  viscous  boundary 
layer.  The  two  flows  were  to  be  calculated  simultaneously.  It  was  thought  that 
this  was  a  more  feasible  task  than  trying  to  solve  the  Reynolds  averaged  Navier- 
Stokes  equations  directly.  Because  of  the  important  rotational  effects  introduced 
by  the  blade  row  when  large  turning  is  involved,  the  Euler  equations  were  used 
to  describe  the  inviscid  parts  of  the  flowfield.  Unfortunately,  it  was  found 
at  the  beginning  of  the  present  program  that  the  state  of  the  art  for  calculating 
numerical  solutions  to  the  Euler  equations  for  the  flow  through  blade  rows  was 
not  sufficient  to  support  the  overall  goals  of  the  program.  In  fact,  two- 
dimensional  and  three-dimensional  Euler  code  solutions  for  cascades  are  only 
recently  emerging  in  the  open  literature.  Furthermore,  the  majority  of  these 
schemes  are  explicit  time  marching  schemes  that  are  extremely  long  running  codes 
that  make  it  impractical  to  incorporate  an  interacting  boundary  layer. 

The  approach  that  was  followed  in  the  present  program  to  calculate  the 

inviscid  flow  was  to  complete  the  development  of  an  implicit  time  marching  Euler 
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code  that  was  initiated  under  AFOSR  sponsorship.  ’  The  two-dimensional  version 
of  this  code  has  now  been  developed  for  cases  of  subsonic  relative  inflow  where 
either  a  Kutta  condition  is  specified  at  the  blade  trailing  edges  or  the  condi¬ 
tions  at  downstream  infinity  are  specified.  The  inclusion  of  a  Kutta  condition 
proved  to  be  particularly  troublesome  and  consumed  so  much  effort  that  inclusion 
of  an  interacting  boundary  layer  was  not  feasible  under  the  present  program. 
However,  now  that  this  problem  is  solved,  the  extension  of  the  code  to  three 
dimensions  is  straightforward  although  a  viable  algorithm  for  incorporating  an 
interactive  boundary  layer  is  yet  to  be  developed. 

The  work  on  analysis  of  the  unsteady  flow  through  a  blade  row  was 

initially  approached  from  two  viewpoints.  The  specific  problem  of  interest  was 

the  response  of  a  blade  row  to  inlet  distortion.  The  first  approach  was  the 

completion  of  the  linearized,  unsteady,  lifting-surface  theory  that  was  initiated 

2 

under  a  previous  program.  Also  considered  was  the  development  of  a  nonlinear 
transonic  theory  that  would  account  for  finite  blade  thickness  and  loading  effects. 
However,  after  some  study  it  was  realized  that  the  time  marching  Euler  code  being 
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developed  for  the  prediction  of  steady  state  performance  could  also  handle  the 
unsteady  distortion  problem  if  the  boundary  conditions  at  infinity  and  the 
spatially  periodic  boundary  conditions  are  modified  to  show  the  appropriate  time 
varying  behavior.  Hence  this  work  was  expended  on  the  Euler  code  development  and 
every  effort  was  made  to  keep  the  resulting  code  time  accurate  (i.e.,  to  make 
the  calculations  follow  the  actual  flow  development  with  time) . 

The  Euler  code  development  is  described  in  Section  II.  The  two- 
dimensional  stability  theory  work  is  described  in  Section  III  and  the  unsteady 
lifting  surface  theory  is  described  in  Section  IV. 


4 


SECTION  II 

EULER  CODE  DEVELOPMENT 


1 .  INTRODUCTION 

The  present  algorithm  that  has  been  developed  to  calculate  the  inviscid 
flow  through  a  blade  row  consists  of  two  computer  codes.  The  first  code  generates 
a  boundary-fitting  coordinate  system  that  maps  the  physical  blade  row  and  flow 
under  consideration  into  a  box- like  computational  domain.  This  greatly  faci¬ 
litates  the  application  of  blade  surface  boundary  conditions  and  the  specifica¬ 
tion  (or  calculation)  of  flow  conditions  at  infinity  upstream  and  downstream  of 
the  blade  row.  The  basic  mapping  method  used  is  a  modification  of  one  due  to 
Ives  and  Liutermoza.  The  modifications  concern  the  final  destination  of  the 

blade  surface  and  the  points  at  infinity  which  were  assigned  similarly  to 
7  8 

Dulikravich.  ’  The  mapping  code  writes  the  metrics  on  a  tape  which  is  input 
into  the  second  code.  The  second  code  integrates  the  transformed  Euler  equations 
in  the  computational  domain  and  applies  the  appropriate  boundary  conditions. 

The  basic  integration  method  is  an  implicit  time  marching  scheme  that  uses 
approximate  factorization.  The  steady  state  solution  is  obtained  as  the  long 
time  limit  of  an  unsteady  problem  with  steady  boundary  conditions.  This  inte¬ 
gration  scheme  was  developed  for  external  aerodynamic  problems  and  is  most 
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clearly  described  in  a  series  of  papers  by  Beam  and  Warming.  ’  This  integra¬ 
tion  scheme  was  selected,  originally,  because  it  was  judged  to  be  the  most 
promising  in  terms  of  computational  efficiency  and  running  time  of  any  of  the 
methods  then  available  and  it  is  still  thought  to  be  a  leading  candidate  from 
this  viewpoint.  The  explicit  time  marching  codes  require  much  longer  running 
times  because  of  the  CFL  stability  requirements  on  time  step  size. 

The  development  of  both  codes  used  in  the  present  program  was  initiated 
under  an  AFOSR  sponsored  program  and  the  status  of  the  codes  at  the  end  of  that 
program  is  reported  in  Refs.  4  and  5.  At  that  time,  the  integration  code 
was  not  converging  and  the  reason  suspected  was  the  singular  behavior  of  the 
metrics  of  the  transformation  at  the  blade  trailing  edge.  This  problem  was 
solved  and  made  little  difference  in  the  results.  The  critical  points  turned 
out  to  be  the  analytical  method  that  was  being  used  to  calculate  the  metrics 
in  general,  the  time  differencing  scheme  being  used,  and  the  method  of  evaluating 
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quantities  on  the  blade  surfaces.  In  addition,  after  the  code  was  revised  to 
include  these  points  and  produced  converged  results,  it  was  found  that  a  revised 
treatment  of  the  numerical  viscosity  term  eliminated  undesirable  spatial  oscilla¬ 
tions  that  appeared  in  the  initial  calculation.  All  of  these  points  have  been 
cleared  up  satisfactorily  and  the  integration  code  works  over  a  range  of 
subsonic  relative  inlet  conditions  with  either  a  Kutta  condition  specified  or 
the  conditions  at  downstream  infinity  specified.  Global  conservation  of  mass 
flow  is  not  a  problem  with  the  present  code  in  contrast  to  other  published 
results  of  Euler  codes. 

The  overall  algorithm  developed  here  is  time  accurate  except  for  the 
treatment  of  the  inlet  and  outlet  conditions.  This  feature  was  deliberately 
maintained  so  that  the  code  would  only  require  slight  modification  in  order  to 
consider  the  response  of  the  blade  row  to  unsteady  inlet  distortion.  Thus, 
several  techniques  that  have  been  used  on  time  marching  schemes  to  speed  the 
approach  to  the  steady  state  limit  have  been  avoided  in  order  to  maintain  the 
time  accuracy. 

The  basic  analysis  upon  which  the  present  codes  are  based  is  given  in 
Refs.  4  and  5.  In  order  to  fully  explain  the  critical  points  mentioned 
above  and  to  give  a  coherent  description  of  the  existing  codes,  it  will  be 
necessary  to  duplicate  some  of  the  material  contained  in  the  earlier  references. 
Also,  because  the  coordinate  mapping  used  influences  the  manner  in  which  the 
boundary  conditions  are  applied,  it  is  not  possible  to  separate  totally  the 
discussion  of  the  mapping  from  the  discussion  of  the  integration  scheme.  Thus 
the  basic  part  of  the  discussion  will  be  dealing  with  the  Euler  equations  written 
in  a  general  curvilinear  coordiate  system  with  a  subsection  on  the  particular 
coordinate  transformation  that  was  used.  In  addition,  the  basic  analysis  was 
originally  performed  for  three-dimensional  flow  and  this  point  of  view  is  con¬ 
tinued  here  even  though  only  the  two  dimensional  analysis  has  been  fully  imple¬ 
mented  into  the  code. 
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is  the  density 
is  static  pressure 

are  respectively  the  axial,  radial  and  circumferential  velo¬ 
cities  in  the  system  fixed  to  the  blades 

is  the  rotation  rate  of  the  blade  row 
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A  general  transformation  of  geometric  coordinates  is  introduced 


where  the  surface  ~§-  Q  will  eventually  be  the  blade  surfaces  and  >7  will  be 
the  radial  variable.  The  Jacobian  of  the  transformation  is  defined  as 
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where  the  subscripts  stand  for  partial  differentiation. 
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Following  Viviand,  it  is  possible  to  redefine  the  dependent  variables 
in  such  a  fashion  that  the  transformed  Euler  equations  retain  the  same  semi¬ 
conservation  form  as  Eq.  (1).  This  is  accomplished  by  defining  the  following 
quantities : 


G  =  "  -SrF  +  ^ J 


H  -  F 

r  Q 


Then  Eq.  (1)  becomes 

c)  U  +  3  £  *  £  F"  -b  c)  G  -  /-/ 

^  t  £  5  ^ 

9 

Next,  following  Beam  and  Warming  ,  the  time  differencing  scheme  is 

t<  >> 

considered.  Let  the  superscript  n  denote  the  time  step;  then  Taylor  series 
expansion  may  be  used  to  show 
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When  «£*  j>  ,  this  is  known  as  trapezoidal  differencing,  when«*-0  it  is  known  as 
Euler  implicit,  when  <x  -  /  it  is  known  as  Euler  explicit.  Clearly  when 
the  greatest  time  accuracy  is  achieved  as  seen  from  the  error  term.  Eq.  (2) 
will  eventually  be  used  to  replace  the  time  derivatives  which  appear  in  the 
right-hand  side  of  Eq.  (3);  but,  some  preliminary  work  is  required  in  order 
to  end  up  with  a  scheme  that  minimizes  the  computations.  First,  use  is  made  of 
the  homogeneity  property  of  the  matrices  appearing  in  Eq.  (2),  that  is,  they 
may  be  expressed  as 

-A  ^ 

£  -  A  U 

__ ^  -*S 
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•'■'S  -A  -A 

H  =  D  U 
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This  property  allows  a  time  linearization,  at  least  to  the  order  of  Eq.  (3),  for 
then  by  Taylor  series  expansion  one  obtains 
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Introducing  the  new  independent  variable  AU  -  U  ~  U  the  combination 
of  Eqs.  (1),  (3)  and  (4)  can  be  put  in  the  so-called  "delta  form" 
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where  I  is  the  identity  matrix  and  differentials  operate  on  all  products 
formed  to  their  right.  That  is,  the  first  differential  expression  appearing  on 
the  left  hand  side  of  Eq.  (5)  should  be  interpreted  as  ^  (_  A  A  U  )  .  Adding 
terms  of  order  Azf'3to  the  left  hand  side  of  Eq.  (5),  the  following  approximate 
factorization  is  achieved: 


(*  C 

(i  *  (  I-<k)a±(  ^  -Dn)'j  A  U“  *  =  AU  **n 

(i  -  (!-*)  At  Cn)AUn  *  AU*n 
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(6b) 


(6c) 


This  system  has  the  same  time  accuracy  as  Eq.  (5) . 
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This  factorization  effectively  decouples  the  original  three-dimensional 
equation  (Eq.  (5))  into  a  sequence  of  three  one-dimensional  equations  (Eq.  (6)) 
which  allows  a  great  reduction  and  simplification  of  the  following  computations. 
When  central  differencing  is  used  for  the  spatial  derivatives,  second  order 
accuracy  in  the  spatial  variables  is  obtained.  However,  it  has  generally  been 
found  necessary  to  add  numerical  viscosity  terms  to  this  system  of  equations 
in  order  to  prevent  the  occurence  of  spatial  oscillations  which  can  destroy  the 
validity  of  the  solution.  The  actual  system  of  equations  solved  is  then 
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Where  €A  is  the  implicit  numerical  viscosity  coefficient  and  is  the 

explicit  numerical  viscosity  coefficient.  In  the  application  of  this  method 
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to  external  aerodynamic  problems,  5  all  of  the  explicit  damping  terms  were 
appended  to  the  first  factorization.  However,  better  results  were  obtained  in 
the  present  investigation  by  associating  the  explicit  damping  term  in  a  given 
direction  with  the  approximate  factorization  in  that  direction. 

Generally  the  numerical  viscosity  coefficients  are  chosen  to  be  of 
order  At  with  £ A  chosen  to  be  twice  €a  after  the  stability  analysis  performed 
by  Desideri.  When  central  difference  expressions  are  used  for  the  spatial 
derivatives,  Eqs .  (7a),  (7b)  and  (7c)  each  represent  a  block  tri-diagonal 
system  of  equations  which  is  solved  using  the  Thomas  algorithm.  The  order 
of  each  system  and  the  number  of  times  it  must  be  solved  depend  upon  the  size 
of  the  computational  grid  used.  The  subcase  of  two  dimensional  flow  is  obtained 
by  ignoring  the  factorization  in  the  preceeding  analysis. 

3.  TRANSFORMATION  AND  METRICS 

As  previously  mentioned,  the  purpose  of  transforming  coordinate  systems 
is  to  map  the  actual  complex,  three-dimensional  geometry  of  the  blade  row  and  flow 
field  into  a  box- like  computational  domain  where  blade  surface  boundary  conditions 
are  easily  applied. 

With  reference  to  the  previous  sketch,  mapping  starts  from  a  cylindrical 
coordinate  system  in  the  physical  space,  (  ^  Q  ,  /  )  where  /"  is  along  the  blade 
span  and  /  'is  in  the  downstream  axial  flow  direction.  First,  the  geometry  of 
the  axisymmetric  flow  passage  is  used  to  define  a  fractional  distance  between 
hub  and  tip.  Let 


r  -  rsrx) 

rrtx)  -  (x) 


where  r„  (x)  is  the  hub  (inner)  radius  of  the  passage  and  ^r(x)  is  the  tip 
(outer)  radius  of  the  passage.  The  intersection  of  the  blade  surfaces  with 
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surfaces  of  constant  yj  defines  a  two-dimensional  cascade  with  X  and  S  as 
coordinates.  This  cascade  is  ultimately  mapped  to  a  rectangle  in  jf'  space 
as  shown  in  Fig.  1  In  this  figure,  K  is  the  index  in  the  Jf  direction  and 
l_  is  the  index  in  the  'f  direction.  The  figure  depicts  the  case  for  a  10  by 
40  grid  which  was  used  throughout  code  development.  The  mapping  is  patterned 
after  the  work  of  Ives  and  Liutermoza6  with  the  exception  of  the  last  step. 

The  mapping  consists  of  a  sequence  of  six  conformal  transformations  followed  by 
a  shearing  of  coordinates  such  that  the  blades  end  up  along  the  bottom  of  the 
computational  domain  and  the  points  at  infinity  end  up  along  the  top  of  the 

7  8 

computational  domain.  (This  is  a  feature  of  the  mappings  used  by  Dulikravich.  ’  ) 
All  steps  except  the  last  one  are  performed  as  reported  in  Ref.  7. 

The  important  features  of  the  transformation  used  here  are  that  the 
cascade  blade  trailing  edges  and  the  points  at  downstream  infinity  (far  down¬ 
stream  of  the  blade  row)  map  to  the  comers  of  the  computational  domain.  The 
metrics  are  singular  at  these  points.  However,  since  these  points  are  located 
at  the  comers  one  avoids  the  need  to  evaluate  the  difference  equations  contain¬ 
ing  the  metrics.  Any  information  required  in  the  scheme  at  these  points  is 
found  by  extrapolation  or  averaging  from  the  neighboring  points. 

Two  different  cascade  geometries  have  been  used  in  the  program  develop¬ 
ment.  The  first  cascade  blade  is  that  given  in  Ref.  17.  It  is  a  relatively 
fat,  blunt  blade  shape  with  high  solidity,  a~  ,  and  is  here  denoted  as  TW-1. 

The  blade  shape  and  the  boundary  conforming  coordinate  mesh  are  shown  in  the 
physical  plane  in  Fig.  2.  Only  one  blade  is  shown  and  the  width  of  the  mesh 
corresponds  to  the  blade  spacing  in  the  cascade.  This  type  of  boundary  conform¬ 
ing  coordinate  system  is  generally  referred  to  as  an  "0"  coordinate  system 
because  the  $  equal  to  constant  lines  close  on  themselves  downstream  of  the 
blade.  This  blade  shape  is  not  a  practical  aerodynamic  shape  but  was  merely 
chosen  in  Ref.  17  as  a  demonstration  case  whose  ordinates  were  easily  defined. 

It  does,  however,  represent  a  severe  test  for  any  numerical  method  because  of 
its  extreme  thickness  and  high  solidity.  Hence  it  was  retained  here  for  initial 
development  purposes.  In  addition,  a  more  conventional  cascade,  the  NACA  65(12)10 
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with  solidity  of  1  was  chosen  for  preliminary  investigation.  The  ordinates  for 
this  blade  shape  are  given  in  Ref.  18.  Transonic  calculations  for  this  cascade 
are  aslo  given  in  Ref.  16.  The  blade  shape  and  boundary  conforming  coordinate 
mesh  for  this  cascade  in  the  physical  plane  are  shown  in  Fig.  3. 

In  Ref.  5,  the  metrics  were  calculated  analytically  by  application  of 
the  chain  rule  to  the  intermediate  transformations.  As  such,  the  metrics  will 
not  satisfy  the  finite  difference  form  of  Eq.  (7)  exactly  when  the  flow  is 
uniform.  Although  the  error  is  theoretically  second  order  or  higher  in  the 
spatial  differences,  it  produces  a  distributed  source-like  effect  which  evidently 
accumulates  errors  to  result  in  divergence  of  the  calculations  after  only  a  few  . 
time  steps.  This  problem  was  solved  by  calculating  the  metrics  numerically  by 
use  of  the  same  spatial  difference  scheme  that  is  used  to  discretize  the  Euler 
equations.  Then  the  metrics  satisfy  the  equation  exactly  and  the  source-like 
effects  are  avoided.  The  metrics  are  calculated  directly  from  the  coordinates 
themselves  by  using  either  central  differences  or  the  appropriate  one-sided 
differences  along  the  boundary  of  the  computational  domain.  The  appropriate 
expressions  used  are: 

~  ( r>\  ~  r3  ^  <£) 

Vr  ~  (  rsSf  -rf9s)  /£-' 

3.  =  (  /£)-< 
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Vr  r  <  ~ 

=  (0j **>,  - 

¥ Q  ~  (  r rj  )  /  Q 

V*  *  -xsry)  /£" 

~  rf)  /  &  ' 


15 


where 


(x  ,  r,e) 


£> 


-I 


/ 

37 


is  the  Jacobian  of  the  transformation. 

The  problem  with  the  metrics  is  essentially  the  same  problem  as  main¬ 
taining  the  proper  free  stream  conditions  in  external  flow  calculations.  This 
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problem  is  discussed  by  Pulliam  and  Steger,  but  they  indicate  that  the  problem 
is  not  critical  for  external  aerodynamic  calculations.  Evidently,  the  boundary 
conditions  on  the  external  aerodynamics  problem  make  it  a  more  forgiving 
algorithm. 

4.  BOUNDARY  AND  INITIAL  CONDITIONS 

With  the  exception  of  the  comer  grid  points  in  the  transformed 
computational  domain,  initial  conditions  at  t-O  are  needed  at  all  grid  points 
and  boundary  conditions  are  needed  at  each  time  step  for  grid  points  along  the 
boundary  of  the  computational  domain.  First  of  all,  the  boundary  conditions 
will  be  discussed  starting  with  conditions  on  the  blade  surfaces. 

The  general  treatment  of  the  boundary  conditions  follows  Refs.  10 
and  11.  While  evaluating  field  points  &U  is  held  at  zero  on  the  boundaries, 
then  L/  is  explicitly  updated  before  proceeding  on  to  the  next  time  step  calcu- 
lation  of  4< J  ,  Along  the  blade,  surface  flow  tangency  is  imposed  by  considering 
the  components  of  the  velocity  vector  in  the  computational  domain.  They  are 
defined  as 


w,  =  J*  w*  *  fr  Wr  ■* 

Ja  wa 

(8  a) 
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Wj.  =  ■+  y)r  Wr  -+ 
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These  definitions  become  a  system  of  equations  suitable  for  determining  w*  ,  wr  , 
and  Wa  on  the  surface  whqn  lA/3  -O  is  imposed  and  extrapolated  values  of  vJ, 
and  Wz  are  used.  The  value  of  surface  pressure  is  found  from  the  following 
expression  derived  in  Ref.  4. 
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As  explained  in  Ref.  4,  this  equation  is  derived  from  the  momentum 
equations  and  the  continuity  equation  without  benefit  of  the  energy  equation. 

The  equation  has  the  desirable  quality  that  it  does  not  contain  any  time  deriva¬ 
tives.  After  all  of  the  quantities  on  the  right-hand  side  of  the  equation  are 
determined  either  by  the  solution  of  Eqs.  (8)  or  extrapolation,  Eq.  (9)  may  be 
integrated  implicitly  to  determine  surface  pressures.  (In  Ref.  4,  this  equation 
was  integrated  explicitly.  However,  it  was  determined  that  this  was  introducing 
an  instability  into  the  algorithm.)  Then  sufficient  information  is  known  to 
calculate  U  on  the  blade  surfaces.  For  the  subcase  of  two  dimensional  flow, 
we  have  Wz-0  and  Pj?  derivatives  are  neglected.  This  drastically  simplifies 
the  solution  of  Eqs.  (8)  and  (9).  For  three-dimensional  flows,  an  equation 
analogous  to  Eq.  (9)  is  given  in  Ref.  4,  and  O  is  used  in  Eq.  (8)  along 

with  extrapolated  values  for  and  to  provide  sufficient  relations  to 
determine  U  at  the  hub  and  shroud  boundaries. 
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The  remainder  of  this  discussion  deals  with  the  two-dimensional  case. 

The  lines  K  =  1  and  K  =  40  are  the  same  line  in  the  physical  plane  (see  Figs. 

-A 

2  and  3)  and  U  is  determined  on  these  lines  by  averaging  quantities  at  fixed 
L  values  from  the  lines  K  =  2  and  K  =  39  (with  the  exception  of  the  values  at 
L  =  1,  2,  9  and  10  whose  treatment  depends  upon  whether  the  Kutta  condition  is 
applied  or  not) .  The  line  L  =  10  corresponds  to  the  periodic  boundary  in  the 
physical  plane  so  that  these  values  must  be  symmetric  about  the  point  (or 

the  midpoint) .  The  proper  symmetry  is  obtained  by  averaging  from  the  proper 
points  on  the  line  L  =  9.  For  instance,  the  values  at  K  =  3,  L  =  10  are  obtained 
by  averaging  from  the  values  at  K  =  3,  L  =  9  and  K  =  37,  L  =  9  and  assigning  the 
averaged  values  to  the  point  K  =  37,  L  =  10  also.  These  averagings  along  the 
boundary  are  performed  once  for  each  time  step  after  the  interior  field  points 
have  been  calculated. 

The  treatment  of  the  inlet  and  outlet  conditions  that  are  used  in  the 
present  investigation  would  appear  to  constitute  an  overspecification  of  the 
problem  if  the  point  of  view  appearing  in  Refs.  19  and  20  is  adopted.  These 
references  present  an  argument  for  specifying  three  out  of  the  four  dependent 
variables  at  the  inlet  and  one  out  of  the  four  dependent  variables  at  the  outlet. 
These  arguments  are  based  upon  consideration  of  the  origin  of  the  characteristic 
lines  that  pass  through  a  given  point  on  the  upstream  or  downstream  calculation 
boundary.  The  unknown  quantities  at  each  boundary  are  then  to  be  solved  for  by 
using  the  compatibility  relation  along  the  characteristics  that  pass  through 
each  grid  point  on  the  boundary.  It  would  appear  that  the  definitive  analysis 
of  the  proper  boundary  conditions  to  apply  would  involve  applying  a  radiation 

condition  at  each  boundary.  Although  preliminary  analysis  along  these  lines 

21  27 

have  appeared,  ’  “  this  approach  has  not  been  fully  implemented.  Moreover, 
several  methods  of  specifying  inlet  and  outlet  conditions  appear  to  produce 
converged  results.  With  this  situation  prevailing,  it  was  decided  to  investigate 
the  simplest  conditions  first  and  to  choose  the  easiest,  most  computationally 
economical  boundary  conditions  to  apply.  In  this  vein,  the  conditions  that 
apply  in  the  steady  state  limit  have  proved  satisfactory  and  are  described  in 
the  following. 
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In  all  cases,  the  inlet  conditions  at  -  oo  are  completely  assigned.  This 
corresponds  to  assigning  all  the  components  of  U  at  L  =  10,  K  =20  and  21  (a 
common  value  for  all  time  steps) .  In  cases  where  the  outlet  conditions  are 
fixed,  the  outlet  conditions  are  calculated  from  an  assigned  static  pressure  ratio 
and  use  of  the  isentropic  relationship,  the  steady  state  continuity  equation, 
and  the  specified  inlet  conditions.  These  calculated  outlet  values  are  applied 
at  the  points  K  =  1,  L  =  9  and  K  =  2,  L  =  10  and  also  at  the  corresponding  points 
K  =  40,  L  =  9  and  K  =  39,  L  =  10.  In  cases  when  the  outlet  condition  is  fixed, 
no  flow  quantities  are  calculated  at  the  trailing  edges  of  the  blades,  i.e., 
at  L  =  1,  K  =  1  and  40. 

In  cases  when  a  Kutta  condition  is  applied,  the  conditions  at  the  trail¬ 
ing  edge  are  calculated  by  averaging  the  values  at  neighboring  points  on  the 
blade  surface  and  one  step  above  (i.e.,  L  =  2  neighbors).  This  value  is  also 
assigned  at  the  points  L  =  2,  K  =  1  and  40.  The  outlet  condtions  at  K  =  2, 

L  =  10  are  found  by  requiring  no  gradient  in  the  streamwise  direction  to  exist 
at  this  point.  These  values  are  also  assigned  at  K  =  1 ,  L  =  9  and  the  correspond¬ 
ing  points  K  =  39,  L  =10  and  K  =  40,  L  =  9.  In  addition,  the  outlet  mass  flow 
is  preserved  at  each  time  step  by  assigning  y?vvr  the  value  it  has  at  upstream 
infinity. 

The  initial  conditions  used  correspond  to  those  of  uniform  flow  at  the 
inlet  conditions.  The  boundary  conditions  imposed  on  the  blade  surface  also 
correspond  initially  to  the  same  uniform  flow  and  are  gradually  switched  over 
to  the  correct  boundary  condition  using  an  exponential  function  to  smoothly 
transfer  conditions.  Generally  about  50  to  100  time  steps  are  required  before 
the  blade  surface  boundary  conditions  are  fully  applied. 

5.  NUMERICAL  EXPERIENCE  WITH  CODES 

The  initial  modifications  made  on  the  codes  described  in  Refs.  4  and  5 
were  concerned  with  the  metrics,  the  transformation  of  the  points  at  infinity 
and  the  implicit  calculation  of  blade  surface  boundary  conditions.  These  modi¬ 
fications  had  the  effect  of  greatly  increasing  the  number  of  time  steps  before 
divergence  occurred,  but  divergence  still  occurred! 
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Then  the  effect  of  the  time  differencing  scheme  used  was  investigated. 

In  the  original  integration  code,  the  trapezoidal  differencing  was  used  (  QC=  '/£  ) 
because  of  its  greater  time  accuracy.  However,  it  was  found  that  converged 
results  could  be  obtained  at  the  lower  inlet  Mach  numbers  with  Euler  implicit 
differencing  («.  =  Q).  The  results  from  these  first  converged  calculations  are 
shown  in  Figs.  4  through  7.  In  Fig.  4,  the  normalized  surface  pressure  is 
plotted  for  the  TW-1  cascade  for  several  elapsed  times  after  initiation  of  the 
calculations.  The  surface  pressures  are  plotted  as  a  function  of  K,  the  jf 
index,  in  the  computational  plane.  Inspection  of  Fig.  2  shows  that  A  5  is  nearly 
constant  in  arc  length  along  the  blade  surface  so  that  Fig.  4  is  nearly  equiva¬ 
lent  to  unwrapping  the  blade  surface  shape  and  plotting  surface  pressure  against 
arc  length  beginning  at  the  trailing  edge  on  the  upper  surface.  The  surface 
pressure  is  normalized  by  the  product  of  free  stream  density  and  square  of  the 
free  stream  axial  velocity.  The  calculation  starts  with  blade  boundary 
conditions  that  correspond  to  uniform  flow  in  the  physical  plane  [as  if  the 
blade  were  absent) .  The  correct  flow  tangency  boundary  conditions  are  then 
gradually  approached  in  an  exponential  fashion  being  fully  applied  after  the 
first  100  time  steps.  The  time  steps  are  counted  by  the  index  IT  with  a  non- 
dimensional  time  increment  between  steps  of  At  =  .005.  The  real  time  is 
normalized  by  U/c  where  (J  is  the  undisturbed  velocity  in  the  axial 
direction  far  upstream  of  the  blade  row  and  C  is  the  blade  chord.  A  disturb¬ 
ance  is  then  convected  along  one  chord  length  at  free  stream  conditions  during 
a  nondimensional  time  of  unity.  For  the  TW-1  cascade,  there  is  still  a  very  small 
time  dependence  for  the  surface  pressure  at  the  end  points  after  1800  steps, 
but  all  of  the  other  points  have  settled  out. 

It  can  be  seen  that  the  basic  character  of  the  solution  has  been 
established  after  300  time  steps.  The  vast  majority  of  calculation  time  is 
spend  merely  in  refining  the  solution.  For  the  case  shown,  where  the  relative 
inlet  Mach  number  =  0.5,  there  is  evidently  a  shock  on  the  upper  surface, 
being  smeared  over  a  couple  of  mesh  widths. 

The  surface  pressure  results  for  the  NACA  65(12)10  cascade  are  shown 
in  Figs.  6  and  7.  The  relative  inlet  Mach  number  is  0.5  and  the  relative  inlet 
angle  is  45°.  Under  these  conditions,  the  flow  is  subcritical  everywhere  and 
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there  are  no  shocks.  This  probably  results  in  the  quicker  convergence  to  a 
steady  state  solution  compared  to  the  TW-1  cascade.  Attempts  to  obtain  converged 
results  (with  this  early  version  of  the  integration  code)  for  this  blade  row  at 
an  inlet  Mach  of  .76  were  unsuccessful.  This  problem  was  initially  thought  to 
be  due  to  possible  overspecification  of  the  inlet  and  outlet  conditions. 

Further  analysis  of  the  numerical  results  reveals  a  slight  loss  in 

rothalpy , 

I  =  —  -2  *  A./ -rw] 
r-i  />  i  1  '  ' 

at  the  surface  points  near  the  leading  edge.  This  is  worse  for  the  TW-1  blade 

row.  However,  the  biggest  failings  of  the  calculation  are  spatial  oscillations 
that  occur  near  the  trailing  edge  on  the  upper  surface.  These  oscillations  were 
removed  by  revising  the  numerical  viscosity  treatment  as  previously  explained. 

Results  for  the  NACA  65(12)10  cascade  using  the  revised  numerical  vis¬ 
cosity  formulation  are  shown  in  Fig.  8  and  the  results  are  compared  to  those 
obtained  with  the  original  numerical  viscosity  formulation.  It  can  be  seen 
that  the  oscillations  at  the  trailing  edge  were  affecting  the  whole  solution. 

In  addition  to  smoothing  the  results  near  the  trailing  edge,  it  was  found  that 
converged  solutions  could  be  obtained  at  higher  inlet  Mach  numbers  over  a  range 
of  pressure  ratio.  Examples  of  these  calculations  are  shown  in  Figs.  9  through 

11.  These  cases  were  run  to  compare  with  the  results  presented  in  Ref.  19. 

However,  an  exact  comparison  cannot  be  made  because  the  calculation  of  Reference 
19  contains  the  effects  of  a  boundary  layer.  However,  the  present  calculations 
show  similar  trends  to  those  of  Ref.  19.  Moreover,  it  may  be  seen  from  Figs. 

9,  10  and  11  that  the  surface  pressure  distribution  is  extremely  sensitive  to 
the  prescribed  pressure  ratio,  especially  at  the  lower  pressure  ratio  where 
there  may  be  a  shock  on  the  lower  surface. 

6.  KUTTA  CONDITION 

A  number  of  schemes  were  investigated  in  order  to  incorporate  the 
Kutta  condition  in  the  present  algorithm.  These  included  attempts  to  integrate 
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the  Euler  equations  along  the  blade  surface  boundaries  and/or  along  the  K  =  1, 

40  boundaries  in  place  of  the  previously  described  averaging  and  extrapolation 
techniques.  Various  end  conditions  at  downstream  infinity  were  also  investigated. 
Generally  the  integration  schemes  seemed  to  go  awry  because  there  turn  out  to  be 
conflicting  requirements  on  the  metrics  in  order  to  avoid  producing  source- like 
effects  at  the  points  adjacent  to  the  comers  of  the  computational  domain.  This 
dilemma  arises  because  different  one-sided  differencing  techniques  than  those 
used  on  the  field  points  must  be  used  on  the  boundaries  in  order  to  avoid  the 
use  of  data  at  the  comer  points  where  the  metrics  are  theoretically  infinite. 

It  was  found,  however,  that  averaging  along  the  K  =  1  and  40  boundaries 
produced  satisfactory  results  when  the  outlet  mass  flow  was  prescribed.  In 
conjunction  with  these  procedures,  the  conditions  at  K  =  2,  L  =  10  could  either 
be  found  by  extrapolation  or  requiring  no  gradients  in  the  streamwise  direction. 
This  later  condition  was  the  one  selected  for  the  final  code  configuration. 
Calculations  resulting  from  this  version  of  the  Kutta  condition  are  shown  in 
Figs.  12  and  13  for  inlet  Mach  numbers  of  .5  and  .76  respectively.  Also,  as  a 
consistency  check  on  the  code,  the  outlet  conditions  obtained  from  the  case 
when  the  Kutta  condition  was  applied  were  used  as  inputs  for  the  case  when 
downstream  conditions  are  specified.  These  results  are  also  shown  in  Figs.  12 
and  13.  For  an  inlet  Mach  number  of  .5,  the  correlation  of  the  two  calculations 
is  good.  The  only  significant  differences  appear  in  the  vicinity  of  the  trail¬ 
ing  edge  and  as  would  be  expected  the  application  of  a  Kutta  condition  smooths 
the  results  in  this  region.  The  correlation  of  the  two  calculations  at  an  inlet 
Mach  number  of  .76  is  also  good.  In  this  case,  there  is  a  shock  on  the  upper 
surfaces  of  the  blades.  This  shock  is  evidently  very  sensitive  to  any  differences 
in  the  flow  at  the  trailing  edge  as  the  correlation  of  the  two  calculations 
is  slightly  off  in  the  shock  vicinity. 

7.  CONCLUDING  REMARKS 

In  summary,  the  two-dimensional  version  of  the  Euler  code  developed 
works  well  over  a  range  of  subsonic  inlet  conditions  including  cases  where 
shock  waves  form  in  the  blade  passages.  The  code  works  with  either  a  Kutta 
condition  or  conditions  at  downstream  infinity  specified.  The  integration  code 
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requires  approximately  2.5  x  10  seconds  per  grid  point  per  time  step  on  a 
Cyber  750  computing  machine  and  requires  160,000  octal  words  of  central  memory. 

For  the  grid  used  which  was  10  by  40  approximately,  800  to  1200  time  steps  were 
required  to  reach  the  steady  state.  The  cases  with  shock  waves  required  the 
longer  times.  The  size  of  the  time  steps  used  correspond  to  Courant  numbers  of 
5  to  10.  The  upper  limit  for  convergence  was  not  established;  however,  the 
results  for  an  inlet  Mach  number  of  .76  appeared  to  suffer  excessive  rothalpy 
loss  near  the  leading  edge  with  the  larger  time  step  sizes.  The  results  for  an 
inlet  Mach  number  of  .5  did  not  show  such  effects  so  that  computing  times  for 
these  cases  can  probably  be  reduced  substantially. 

Although  the  basic  integration  methods  used  for  the  Euler  euqations 
have  been  well  publicized  in  the  external  aerodynamics  literature,  several  modi¬ 
fications  were  found  to  be  crucial  in  order  to  apply  the  method  to,  internal  flows. 
The  calculation  of  the  metrics  of  the  transformation  proves  to  be  extremely 
important  and  a  revision  of  the  numerical  viscosity  treatment  enlarges  and 
enchances  the  domain  where  converged  solutions  can  be  obtained.  In  particular, 
it  was  found  that  the  metrics  must  be  discretized  in  the  same  spatial  fashion  as, 
the  governing  partial  differential  equation  in  order  to  avoid  introducing  source¬ 
like  terms  which  would  quickly  destroy  the  solution. 
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SECTION  III 

TWO  DIMENSIONAL  STABILITY  THEORY 


1.  INTRODUCTION 

A  two  dimensional  stability  theory  to  predict  the  rotating  stall 
inception  boundary  for  incompressible  flow  through  an  isolated  blade  row  and 
stage  (two  blade  rows)  was  presented  in  Ref.  24  and  25.  Correlations  of  the 
theory  with  experiment  on  a  number  of  blade  rows  and  stages  show  that  if  the 
steady  state  loss  and  turning  performance  for  each  blade  row  of  interest  is 
known,  the  inception  boundary  for  rotating  stall  can  be  predicted  quite  accurately. 
Basically  the  theory  shows  that  geometric  parameters  such  as  solidity,  blade 
chord,  stagger  angle,  etc.  influence  the  inception  boundary  and  propogation  speed 
via  the  fashion  that  these  parameters  influence  the  steady  state  performance  of 
the  blade  row  or  stage.  In  the  case  of  a  stage,  the  blade  row  spacing  determines 
the  number  of  stall  cells.  Although  the  blade  row  performance  data  required  as 
input  by  the  theory  is  often  not  available,  especially  for  compressible  flow  * 
conditions,  the  extension  of  the  theory  to  subsonic  compressible  flow  has  been 
accomplished  so  that  qualitative  trends  may  be  obtained  by  using  the  incompressible 
blade  row  performance. 

The  extension  of  the  theory  to  a  single  blade  row  in  compressible  flow  * 
was  given  in  Ref.  26  and  revised  in  Ref.  3.  The  revised  theory  was  extended  to 
two  blade  rows  and  numerical  results  were  obtained  under  the  present  program. 

The  results  of  these  analyses  will  be  given  in  the  following  sections. 

2.  TWO  BLADE  ROW  ANALYSIS 

The  same  general  type  of  analysis  as  used  in  the  single  blade  row  case^ 
has  been  applied  to  the  two  blade  row  problem.  The  blade  rows  are  modeled  by 
finite-thickness  actuator  sheets  as  shown  in  Fig.  14.  There  are  now  three  flow 
regions  designated  by  i  =  1,  2  and  3,  and  the  flow  quantities  are  broken  into 
mean  time- independent  and  time-dependent  parts  as  follows: 
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Axial  velocity:  (J.  (  X  ,  ,  t )  =  Uj_  + 

Tangential  velocity:  1a p(Xt^,t)  =  Wj  +  uJt'  (x  ,  y ,  t ) 

Pressure:  P.  Lx  ,  Lj  ,  t )  -  PL  +  -fs'  (x,ij,t) 

Density:  t")  =  P -L  +  (X  ,  ,  t) 

The  mean  quantities  are  assumed  constant  in  each  flow  region.  The 
absolute  swirl  is  denoted  by  S(  and  is  equal  to  \^J. /U’  •  The  swirl  relative 
to  the  first  blade  row  is  denoted  by  and  is  equal  to  (  W/ -  V4^  )/  UL  ■ 

The  swirl  relative  to  the  second  blade  row  is  denoted  by  ^  and  is  equal  to 
(  W;  -  W5  )/  .  The  disturbance  quantities  are  assumed  to  have  a  time  and 

dependency  of  the  form  expj(  ct  +  >1  ^>/r  ).  These  are  "  n  "  celled  waves 
propagating  around  the  annulus.  Here  C  =  C^jCr  and  the  flow  is  neutrally 
stable  when  Cx  -  O  ,  unstable  if  Cr  <  O  and  stable  if  Cz>0  ■  The  time- 
dependent  quantities  are  assumed  to  be  much  smaller  than  their  corresponding 
steady  parts.  Under  these  assumptions,  the  equations  of  motion  may  be  linearized 
and  the  disturbance  quantities  have  the  following  forms: 
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where 


a.  - 


speed  of  sound  in  flow  region  c 
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Here  3/  ,  and  are  constants  to  be  determined  by  the 

matching  conditions.  The  terms  associated  with  represent  vorticity 

waves,  the  term  associated  with  £ i  represents  an  entropy  wave  and  the 
remaining  terms  represent  potential  (irrotational)  disturbances.  Requiring 
irrotational ,  isentropic  flow  upstream  gives.  D^=E.-0  .  Requiring  bounded 
flow  conditions  upstream  and  downstream,  results  in  P ,  -6^-0  .  These 

considerations  reduce  the  number  of  unknown  constants  to  eight.  These  remain¬ 
ing  unknowns  are  determined  by  applying  four  matching  conditions  across  each 
actuator.  These  matching  conditions  are: 


Conservation  of  Mass  Flow 


Ri+i  “-L  +  i  +  Pui  UL  +  1  -Pi  Ui  +  ***&,<£>  i  +  ,  +  Pi  )  =  0  (11) 
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Blade  Row  Turning  Relation 
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Vorticity  Compatibility  (derived  from  the  momentum  equation) 
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Conservation  of  Energy 
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In  these  relations,  (V  -  )  is  the  loss  coefficient  for  blade  row  4  and 

L  t* 

GR  sd  ■)  is  the  turning  relation  for  blade  row  i  .  A  prime  on  a  quantity 
indicates  its  derivative  with  respect  to  its  argument.  These  relations  are 
all  applied  in  a  blade  fixed  coordinate  system.  A  dot  over  a  quantity  stands 
for  its  time  derivative  and  is  the  vorticity  in  region  i  .  In  order  to 
apply  these  relations,  c  is  set  equal  to  the  actuator  number  and  Equations  la¬ 
id  are  used  to  calculate  the  required  quantities.  At  actuator  1,  flow  quantities 
with  a  subscript  1  are  evaluated  at  X-O,  y  -  and  quantities  with  a  subscript  2 
are  evaluated  at  X  =  J.1  c<x>  S1  ,  ij  -  <j  +A1Mr*6'i.  At  actuator  2,  flow  quantities  with 
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a  subscript  2  are  evaluated  at  and  quantities  with  a  subscript  3  are 

evaluated  at  TC-  T  +  hj  Sz  ,  y-y  -hJ'zAt^  Sz  •  The  result  of  performing  all 
these  operations  is  a  matrix  equation: 
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where  the  fi  matrix  is  8  x  8.  This  equation  has  non-trivial  solution  only  when 
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This  is  the  characteristic  equation  which  determines  the  value  of  C  for  a  given 
flow  configuration  and  blade  row.  After  factoring  out  common  exponential  factors, 
the  expression  for  the  R  are  as  follows: 
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Inspection  of  the  eighth  column  of  the  matrix  shows  that  one  root  of 


Equation  (15)  is  O  or 


C  = 


-  +  i 


2  U3 


J.  6 

2  c* 


(16) 


This  is  a  stable  wave  fixed  to  the  second  blade  row.  This  root  is  sub¬ 
sequently  factored  out  of  and  Rgg  .  Inspection  of  the  n ij  terms  shows 

that  Eq.  (15)  is  transcendental  in  C  and  therefore  must  be  solved  numerically. 

The  analysis  for  a  single  blade  row  may  be  considered  as  a  subcase  of 
the  above  analysis.  For  the  single  blade  row  the  characteristic  determinant, 
/A!  >  is  four  by  four  and  its  elements  may  be  obtained  from  those  presented  by 
ignoring  the  Acj  for  7  4-  and  then  factoring  out  common  exponential  terms. 
Also,  for  the  single  blade  row  case  one  root  can  be  determined  analytically  as 

Ut  =  "'2/5ec  6‘ 

This  root  represents  a  wave  fixed  to  the  blade  row  and  is  always  stable.  The 
remaining  stability  determinant  is  also  transcendental  in  C  and  must  be  solved 
numerically. 


3. 


RESULTS 


Computer  codes  were  written  to  solve  both  the  single  blade  row  and  two 
blade  row  characteristic  equations.  Both  programs  are  quite  similar;  the 
characteristic  equation  is  solved  iteratively  using  a  Newton-Raphson  scheme.  The 
initial  guess  required  at  the  first  set  of  inlet  conditions  is  obtained  from 
the  incompressible  result  for  the  first  (or  only)  blade  row.  After  the  solution 
has  been  obtained  at  this  inlet  condition,  the  inlet  angle  is  then  increased  by 
a  small  increment  and  the  first  solution  is  used  as  the  initial  guess.  The 
inlet  angle  is  then  incremented  until  the  solution  is  obtained  over  the  range 
required.  The  single  blade  row  program  has  been  applied  to  two  blade  sets.  Rotor 
No.  1  and  Stator  Set  No.  4  of  Ref.  25.  Both  blade  rows  were  run  with  inlet  condi¬ 
tions  corresponding  to  sea  level  standard  conditions  and  an  inlet  axial  velocity  of 
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60  ft/sec.  These  conditions  approximate  closely  the  actual  tested  conditions. 

The  upstream  axial  Mach  number  is  0.057,  which  is  essentially  incompressible 
flow.  The  resulting  stability  curves  and  propagation  velocities  agree  well  with 
the  predictions  for  incompressible  flow.  Stator  Set  No.  4  has  also  been  analyzed 
with  upstream  axial  flow  Mach  numbers  of  0.2  and  0.4  using  the  turning  and  loss 
performance  from  incompressible  flow.  These  calculations  are  shown  in  Figs.  15 
and  16  where  they  are  compared  with  the  results  of  our  previous  compressible 
flow  stability  theory  from  Ref.  26  (Figs.  37  and  38).  The  solid  lines  in  both 
figures  are  the  new  calculations  and  the  symbols  are  the  calculations  from 
Ref.  26.  The  difference  between  the  analysis  presented  in  Ref.  26  and  the 
present  analysis  lies  in  the  energy  matching  condition  used.  The  present 
analysis  eliminates  consideration  of  the  viscous  dissipation  function  within 
the  blade  row.  This  function  was  evidently  numerically  very  small  because  it 
scales  with  the  viscosity  and  the  condition  under  consideration  corresponds  to 
relatively  high  Reynolds  number.  The  smallness  of  the  dissipation  function 
evidently  is  responsible  for  the  close  numerical  agreement  in  stability  bounda¬ 
ries  for  the  two  forms  of  the  energy  matching  condition. 

The  calculations  from  the  present  single  blade  row  program  for  Rotor 
No.  1  are  shown  in  Figs.  17  and  18.  The  incompressible  loss  and  turning  perform¬ 
ance  were  used  for  all  the  Mach  numbers  shown.  These  calculations  then  give  the 
influence  of  compressibility  in  the  flow  outside  the  blade  row  on  the  rotating 
stall  boundary;  the  effect  of  compressibility  on  blade  row  turning  and  loss 
performance  is  not  included.  For  both  the  rotor  and  the  stator  these  compressi¬ 
bility  effects  are  seen  to  reduce  the  stability  of  the  flow  with  only  a  small 
change  in  the  neutral  stability  boundary.  Likewise,  the  propagation  velocities 
for  the  stator  set  are  uniformly  reduced.  The  propagation  velocities  for  the 
rotor  show  a  mixed  effect. 

The  two-blade-row  program  was  run  to  correlate  with  the  high  hub-to-tip 
ratio  rotor-stator  tests  reported  in  Ref.  3.  The  rotor  had  a  stagger  angle  of 
-40.0  degrees  and  the  cases  for  stator  stagger  angles  of  28.2  degrees  and  37.2 
degrees  were  considered.  The  rotor  had  a  chord  of  0.12057  feet  and  the  stator 
chord  was  0.10833  feet.  The  axial  velocity  for  the  tests  was  59.46  feet  per 
second  (  /Vj  =  .  05 3) . 
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As  a  check  out  of  this  program,  the  input  data  were  adjusted  such  that 
one  of  the  two  blade  rows  did  not  influence  the  flow.  That  is,  one  of  the  blade, 
rows  did  not  turn , the  flow  or  introduce  any  losses;  this  blade  row  will  be 
referred  to  as  the  inactive  blade  row.  Although  this  blade  row  did  not  influence 
the  flow,  the  flow  is  still  divided  into  three  regions  and  matching  conditions 
are  still  applied  at  both  blade  row  locations  so  that  the  full  program  is 
exercised  in  calculating  the  stability  of  this  case.  The  turning  and  loss 
performance  of  the  rotor  (in  the  presence  of  the  stator)  was  used  for  the  blade 
row  which  influenced  the  flow.  The  results  of  these  calculations  for  a  stator 
stagger  angle  of  37.2  degrees  are  shown  in  Fig.  19.  The  solid  curve  in  this 
figure  is  the  stability  curve  for  the  stage  as  presented  in  Fig.  23,  Ref.  26. 

The  dashed  curve  is  the  stability  curve  for  the  rotor  alone  but  with  rotor  turn¬ 
ing  and  loss  performance  as  measured  in  the  presence  of  the  stator.  Both  these 
curves  are  calculated  using  the  incompressible  flow  stability  theory  from  Ref. 

25.  The  symbols  represent  values  calculated  from  the  present  compressible 
stability  program  with  one  or  the  other  of  the  blade  rows  inactive  as  noted 
and  with  an  axial  velocity  of  59.46  feet  per  second.  The  agreement  with  the 
rotor-alone  calculations  for  incompressible  flow  is  quite  good  and  validates 
the  code. 


The 'two-blade  row  program  was  then  used  to  investigate  the  effects  of 
inlet  axial  Mach  number  and  blade  row  spacing  assuming  that  these  parameters  do 
not  inlfuence  the  blade  row  turning  and  loss  performance  of  either  rotor  or 
stator  (in  the  absence  of  the  appropriate  blade  row  data  under  compressible 
flow  conditions) . 

The  influence  of  axial  Mach  number  on  the  stability  characteristics  of 
the  rotor-stator  stage  with  stator  stagger  angle  of  28.2  degrees  is  shown  in 
Fig.  20.  There  the  results  are  for  one  stall  cell  (n  r  j)  which  has  the  lowest 
stability.  The  neutral  stability  point,  where  the  damping  factor  goes  to  zero, 
is  slightly  decreased  in  magnitude  by  increasing  axial  Mach  number.  However, 
the  stability  in  the  stable  range  of  inlet  swirls  is  drastically  affected.  The 
influence  of  compressibility  seems  to  drive  the  stage  stability  curve  closer 
to  the  rotor  alone  curve  near  the  neutral  stability  point.  Similar  calculations 
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are  shown  in  Fig.  21  for  the  rotor-stator  with  stator  stagger  angle  of  37.2 
degrees.  For  this  particular  geometry,  increasing  axial  Mach  number  has  a 
more  pronounced  effect 'over  the  entire  range  of  inlet  swirl.  The  magnitude  of 
the  inlet  swirl  for  the  neutral  stability  point  decreases  markedly  with  axial 
Mach  number,  falling  approximately  18  percent  between  M,  =  0  and  A \-C.3  .  Again, 
the  shape  of  the  neutral  stability  curve  for  the  stage  tends  toward  that  of  the 
stability  curve  for  the  rotor  alone  as  axial  Mach  number  is  increased.  The 
greater  influence  of  compressibility  on  this  case  evidently  occurs  because  the 
rotor  alone  becomes  unstable  at  a  greatly  reduced  inlet  swirl  magnitude  compared 
to  the  stage. 

The  effects  of  blade  row  spacing  for  the  rotor-stator  stage  with  stator 
stagger  angle  of  28.2  degrees  has  been  investigated  at  an  axial  Mach  number  of 
0.4.  (Again,  the  blade  row  turning  and  loss  performance  is  assumed  not  to  vary 
with  Mach  number  and  also  blade  spacing  in  this  case.)  The  results  are  shown 
in  Fig.  22.  In  this  figure,  T  is  the  distance,  in  feet,  between  the  leading 
edge  of  the  rotor  and  leading  edge  of  the  stator.  The  value  of  T  for  the 
configuration  tested  was  7"  =  .12196  resulting  in  a  gap  of  .02952  feet  between 
trailing  edge  of  the  rotor  and  leading  edge  of  the  stator.  This  corresponds  to  a 
ratio  of  gap  to  rotor  chord  of  .244.  Calculations  were  made  with  values  of  7” 

(in  feet)  of  .1823,  .2426  and  1.823  corresponding  to  ratios  of  gap  to  rotor 
chord  of  .744,  1.244  and  14.34.  The  calculation  for  7”  =  .1823  are  not  shown 
as  they  were  only  slightly  less  than  those  for  7""  =  .2426.  It  is  seen  that 
increasing  the  gap  initially  is  slightly  stabilizing;  however,  as  the  gap  is 
made  very  large,  the  results  approach  the  rotor  alone  case. 

The  axial  wavelength  of  the  rotating  stall  wave  has  been  calculated  for 
purposes  of  comparing  the  blade  spacing  to  the  axial  wavelength  of  the  rotating 
stall  wave  at  neutral  stability  conditions.  For  this  case,  there  are  two  axial 
wavelengths;  one  associated  with  the  rotational  part  of  the  disturbance,  and  one 
associated  with  the  irrotational  part.  At  neutral  stability  conditions,  the 
irrotational  portion  has  a  wavelength  of  approximately  17.4  feet  and  the 
rotational  part  has  a  wavelength  of  approximately  3.2  feet.  It  is  seen  then  that 
the  spacings  investigated  in  Fig.  22  are  generally  much  smaller  than  either  of 
the  axial  wavelengths  associated  with  the  stall  wave. 


37 


4. 


CONCLUDING  REMARKS 


A  small  disturbance  stability  analysis  has  been  presented  for  the  two- 
dimensional,  compressible  flow  through  a  rotor-stator  combination.  The  theory 
has  been  implemented  into  a  computer  code  which  predicts  the  stability  boundary 
and  propagation  speed  of  allowable  disturbance  for  given  blade  row  performance 
data.  The  theory  and  numerical  results  reduce  to  the  incompressible  flow  case 
as  inlet  Mach  number  approaches  zero.  The  calculations  show  that  even  for  fixed 
blade  row  performance,  compressibility  has  a  significant  effect  on  the  stability 
boundary.  It  is  generally  found  that  this  effect  is  destabilizing  compared  to 
the  incompressible  flow  results. 
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SECTION  IV 


THREE-DIMENSIONAL  UNSTEADY  LIFTING  SURFACE  THEORY  FOR  A  ROTOR 


This  portion  of  the  work  is  concerned  with  predicting  the  unsteady 
aerodynamic  response  of  an  annular  blade  row  to  inlet  flow  distortion.  A 
linearized,  but  fully  three-dimensional,  lifting-surface  analysis  was  chosen 
for  this  purpose.  With  small  modifications,  the  analysis  could  be  applied  to 
the  study  of  blade  flutter  as  well. 

1.  FLOW  MODEL 

The  same  geometry  is  assumed  for  the  present  case  as  had  been  used  in 

2 

our  earlier  study  of  the  steady  load  problem.  That  is,  the  blade  row  is  assumed 
to  be  housed  in  an  infinitely  long,  hard-walled  annular  duct  of  constant  hub/tip 
ratio,  Vr  ,  as  shown  in  Fig.  23.  The  flow  is  assumed  to  be  inviscid  and  to 
contain  a  uniform  subsonic  axial  component  at  Mach  number  M  =  U  /  ao  ,  where  &o 
is  the  undisturbed  sound  speed.  Any  inflow  distortions  are  viewed  as  small 
perturbations  about  this  undisturbed  state.  The  blades  rotate  with  a  constant 
angular  velocity,  which  in  this  analysis  is  denoted  _Q  ,  rather  than  u>  as  used 
previously.  In  keeping  with  the  usual  convention,  the  latter  symbol  will  be 
used  for  the  harmonic  time  dependence  below.  Since  the  blade  boundary  condi¬ 
tions  are  more  easily  expressed  in  blade-fixed  coordinates,  we  again  express  the 
governing  equations  in  these  terms.  In  this  frame,  the  time-averaged  undisturbed 
inflow  has  a  velocity  U R  -  [_UZ  +  (n.R)Z]}  and  follows  the  helical  stream  surfaces 
defined  by  7  =  ^  "  fr  z  =  constant.  The  unsteady  flow  is  assumed  to  be  a  small 
perturbation  about  this  flow,  so  the  disturbance  field  will  be  irrotational  and 
isentropic.  The  linearization  again  allows  us  to  apply  the  blade  boundary  condi¬ 
tions  along  the  undisturbed  stream  surfaces,  so  that  to  first  order  in  the 
perturbation  scheme,  blade  thickness  and  camber  do  not  affect  the  unsteady  loads. 
We  also  assume  that  the  undisturbed  relative  Mach  number  is  subsonic  all  along 
the  span. 
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Under  the  above  assumptions,  one  can  derive  a  linear  integral  equation 
relating  the  prescribed  normal  velocity  distribution  on  the  blade  surfaces  to 
the  unknown  loading,  which  appears'  under  the  integral  .  The  steps  involved  have 
already  been  described  in  some  detail  in  Ref.  2,  and  the  final  result  is  quoted 
there  as  Eq.  (63);  because  of  its  length,  it  is  reproduced  here  in  Appendix  B. 
This  equation  then  is  the  starting  point  for  the  work  reported 


2. 


LOADING  EXPANSION 


As  anticipated  in  Ref.  2,  we  have  applied  the  same  techniques  to  invert¬ 
ing  this  integral  equation  as  we  had  used  in  the  steady  problem.  The  first  step 
is  to  assume  that  the  unsteady  pressure  difference  across  the  blades  can  be 

k 

expanded  in  the  following  form: 

p  (£  =  O')  -  p(Z?  ~  0+) 


Ap  (<Ta)Xo)  s 


4  Po  u 
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Here  o~  and  X  are  dimensionless  radial  and  axial  coordinates  defined  by: 
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(18a) 


X  - 
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-1  4  X  $  / 


(18b) 


where  r T  is  the  tip  radius,  and  C A  is  the  axial  projection  of  the  blade  chord, 
which  for  simplicity  is  assumed  constant.  We  will  use  (  <r  ,  X  )  to  refer  to 
an  arbitrary  field  point,  and  (  >  Xa  )  to  refer  to  a  source  location  on 

the  reference  blade  surface,  €,  =  0.  The  variable  <f>t  is  related  to  through 
the  relation  ? 

cpa  =  c^x  Xa  (19) 

Finally,  po  refers  to  the  undisturbed  density. 


It  should  be  remembered  that  all  of  the  dependent  variables  carry  an  implied 


harmonic  time  factor,  e. 


i  (J  t 


which  is  suppressed  for  the  sake  of  brevity. 
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The  above  expansion  is  completely  analogous  to  that  used  in  the 
steady  flow  problem.  It  contains  the  appropriate  square  root  singularity  at 
the  leading  edge,  %a  =  -1,  in  the  i  =  1  term.  Moreover,  each  term  vanishes 
identically  at  the  trailing  edge,  %0  =  +1,  thus  automatically  satisfying  the 
Kutta  condition.  Generally  speaking,  at  each  point  there  will  be  a  phase 
shift  between  the  unsteady  loading  and  the  normal  velocity  which  will  depend 

and  Mach  number.  This  was  not 


on  position,  the  reduced  frequency  cu  -  iOCa- 


2  U 

true  in  the  steady  flow  problem,  where  by  definition  everything  is  "in  phase". 
To  account  for  this,  A-p  must  now  be  allowed  to  be  a  complex  quantity,  having 
both  real  and  imaginary  parts.  This  is  reflected  in  the  unknown  coefficients, 

■  ,  which  are  now  complex: 


fib, 


(20) 


(The  symbol  l  will  be  used  as  both  a  subscript  and  as  V-1  ’  ;  the  appropriate 
meaning  should  always  be  clear  from  the  context  in  which  it  appears.)  From 
Eqs.  (17)  and  (20)  we  see  that  we  now  have  a  total  of  2  x  NI  x  NJ  unknowns. 


A  significant  advantage  of  using  the  expansion  in  Eq.  (17)  is  that 
when  it  is  substituted  into  Eq.  (B-l)  ,  analytical  expressions  can  be  derived 
for  both  the  radial  and  axial  integrations  which  result. 


3-  RADIAL  INTEGRALS 

Most  of  the  radial  integrations  which  result  from  the  above  substitu¬ 
tion  have  the  form: 

IL  ■  f  <r)  da-  (2D 

A 

where  ^  is  an  integer  that  can  take  on  positive  or  negative  values.  Here 
Knk  and  denote  the  radial  eigenvalues  and  eigenfunctions  which  are  an 
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outgrowth  of  the  requirement  that  the  radial  velocity  vanish  at  the  hub  and 
2 

tip  walls.  The  integers  n  and  k  are  the  azimuthal  and  radial  mode  numbers, 
respectively.  This  class  of  integrals  is  identical  to  those  which  arose  in 
the  steady  flow  problem;  as  shown  in  Appendix  B  of  Ref.  2,  they  can  be  evaluated 
in  terms  of  Lommel  functions,  for  which  computer  subroutines  were  already 
available . 


There  are,  however,  a  few  radial  integrations  which  appear  in  the 
wake  terms  of  the  governing  equation  which  deviate  from  the  form  quoted  above. 
These  have  integrands  which  differ  from  that  in  Eq.  (21)  by  factors  of: 

=  lr  +  (<PT<r)*]  a  -  7,2,3  (22) 

where  <j>T  =  Hrr/U.  Such  factors  prevent  using  the  analytical  evaluation 
directly  because  now  the  dependence  on  <T  outside  the  R  k  is  not  as  simple  as 
just  being  raised  to  an  integral  power.  However,  Dafcr)  is  a  simple  function 
which  decreases  monotonically  from  hub  to  tip;  this  suggests  that  it  could 
be  accurately  replaced  with  a  polynomial  fit.  In  experimenting  with  such 
fits,  it  was  found  that  an  expansion  in  inverse  powers  of cT  gave  a  better 
representation  than  one  in  positive  powers  of  CT  .  Hence,  we  set 


Ox  (‘T  ) 


Z  A 


27  =  ) 
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-tJ 


(23) 


Fora  =  3,  <fiT  =  2,  and  AA,  =  4  this  was  found  to  represent  D,  (cr)  with  a 
maximum  error  of  approximately  1.7%  over  the  range  0.4  ^  0~  •S-  1.0.  The 
original  function  was  matched  at  the  points  O'  =  0.4,  0.6,  0.8  and  1.0.  This 
is  probably  an  extreme  case;  for  cl  -  1  or  2  and  lower  values  of  4>T  ,  the  fit 
would  be  even  better.  Using  such  a  procedure,  the  radial  integrals  involving 
Do t(cT  )  can  be  represented  as  a  sum  of  integrals  of  the  form  given  in  Eq.  (21) , 
and  thence  evaluated  analytically. 
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4. 


AXIAL  INTEGRALS 


A  greater  variety  of  axial  integrals  occurs  than  was  the  case  in 
the  steady  flow  problem.  As  a  convenience,  define  the  following  shorthand 
for  the  axial  loading  functions  in  Eq.  (17) : 
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The  integrals  which  arise  from  the  wake  terms  in  the  integral  equation  then 
have  the  form: 

C  4  (X.) 

.SLCx ) 

where  and  indicate  real  and  imaginary  parts,  respectively. 

Before  quoting  the  form  of  the  integrals  in  the  exponential  terms, 
a  few  words  are  in  order  regarding  the  notation.  As  above,  a  capital  C  or  S 
will  be  used  to  refer  to  integrals  involving  the  cosine  or  sine  function, 
respectively.  This  is  followed  by  a  lower  case  s  if  in  addition  the  function 
sgn  (  X  ~%a  )  appears  in  the  integrand.  The  subscripts  n  and  k  again  refer 
to  a  specific  combination  of  azimuthal  and  radial  mode  numbers.  The  super¬ 
script  l  refers  to  one  of  the  axial  loading  functions  defined  by  Eq.  (24)  . 
Finally,  a  second  superscript,  A  or  B,  is  used  to  indicate  whether  the  integral 
is  appropriate  to  a  mode  above  or  below  the  cut-off  condition. 
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The  form  of  the  integrals  above  cut-off  is  then 
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Below  cut-off,  the  integrals  take  on  the  form: 
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The  analytical  evaluation  of  these  axial  integrals  can  be  carried 
out  using  exactly  the  same  method  outlined  in  Appendix  C  of  Ref.  2.  Although 
straightforward,  the  algebraic  manipulations  are  tedious;  the  final  results 
are  given  in  Appendix  A.  It  should  be  noted  that  the  number  of  wake  term 
integrals  is  limited  by  the  fact  that  they  do  not  depend  on  the  mode  numbers. 
Also,  for  a  given  (rt  ,  k)  mode,  one  need  evaluate  only  one  or  the  other  of 
the  sets  of  integrals  defined  in  Eqs .  (26)  and  (27),  since  a  given  mode  is 
either  above  or  below  cut-off. 


5.  GOVERNING  EQUATION 

We  can  now  quote  the  new  form  of  the  governing  equation,  with  the 

■A 

above  definitions  and  substitutions  for  Ap  .  Due  to  its  length,  we  will 
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separate  it  into  its  component  parts,  viz.,  the  wake  terms,  the  axisymmetric 
(  n  -  0)  exponential  terms,  and  the  asymmetric  (  n  #  0)  exponential  terms. 

The  wake  terms  are  given  by 


(R&  jj-  (<r,  x) 
Jm  ~  (<T, 


NX  AJ J 


=  L  L 


L  - 1  Ji-1 


-/3r\T<5i°o  v  * 

■4  tv  <r  [i  +  (cpT  <r)z] 


+  s\x)b 


l C  (X)  by  -  S  Wo,- 
l1  1 


+ 


8  r)T  (<t>T  <r)  £ 


'PO 


4  rc  <r  [i  -t-  (<pT<r )z  ] 


C  (X)  a.'j  +  S  (x)  bi ^ 

L  c\x)  bi :  -  s\x)  a; : 

*  uj 


„  »K  1^~Z  R  (0~) 

(i-o2  z 

k  « x  H°k 


+ 


B  r\. 


4  nvr  [i  ■b(4>r<r)z] 


\C(x)  aLj  +  3L(X)  bij 


MM 


NX 

z 


Rnk  c<r) 


ion  =  -WM  k  a-r 
>  <7 


n 2 


C  {X)b  -  S  C%)  n  ^ 

[(i  +  tT*„  <r*)  (j  -i)2  I  i"fc2  -  (1  +  (^r  <rA  4>T*ncr*)  I*k  ~ 


n  k 


8  4>T  u> 

4  rr  cr 


C  i(X)  4  SL  (X)  b^  • 

l  C‘cx)  ^  -  5  4cx;  a 


MM 


2 


k=r  ^  n  r)k 

n  4  o 


r  2  <ro 


1 


L  n  [i  +  (<pT  cra  )z  ] 

+  44, ,(!*:•"-  d2i } 

3^)  '  nk  ~  T’  m 


n 


.  /  ,  ,  _  1-V47  _  i~^Tl 


■k  f*' 

'  -x>  +  3 
r>  /< 


+  <*,,)  «» 

+  Li  ^  nk 

IT  »>-* 


(28) 


46 


where  the  upper  and  lower  lines  in  the  curly  brackets  represent  the  contribution 
to  the  real  or  imaginary  part  of  /UR  .  On  the  right  side,  B  is  the  blade 


number,  and  the  factors  Rnk  are 

defined  by 
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Recall  • that  the  azimuthal  mode  number,  n  ,  is  related  to  the  summa¬ 
tion  index, m  ,  by  (Ref.  2,  p.  38): 

n  =  m  3  +  p  (30) 


where  p  is  the  shift  in  azimuthal  mode  number  produced  by  the  interblade  phase 
shift,  \xQ  ;  viz. , 

8  ^8  (31) 
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P  = 


The  asymmetry  in  the  values  of  n  allowed  by  Eq.  (30)  prevents  us  from  combining 
the  contributions  from  the  -m  terms  with  those  from  the  +m terms,  as  we  were 
able  to  do  in  the  steady  problem. 


The  first  two  lines  in  Eq.  (28)  are  the  axisymmetric  contribution, 
i.e.,  the  n  -  0  modes.  We  see  that  such  modes  can  arise  only  in  the  special 
case  m=  -p  =  0.  For  this  reason,  these  terms  have  been  written  with  the 
Kronecker  delta  symbol,  6(3o  ,  which  is  unity  when  p  =  0,  but  zero  otherwise. 

2 

Note  that  the  terms  proportional  to  (  ^  -  1)  in  the  second  line 
and  the  first  portion  of  the  fourth  line  vanish  when  =  1.  Since  this  is 
the  only  term  in  Eq.  (17)  representing  uniform  radial  loading,  it  can  be  in¬ 
ferred  that  these  terms  arise  from  trailing  vorticity.  This  component  of  the 
wake  vorticity,  also  present  in  the  steady  problem,  has  its  axis  aligned 
with  the  helical  undisturbed  stream,  and  is  produced  by  spanwise  gradients 
in  the  loading.  There  are  other  terms  in  Eq.  (28)  which  remain  for  ^  =  1, 


47 


but  vanish  when  u3  =  0.  These  terms  can  be  inferred  to  arise  from  the  shed 
vorticity,  whose  axis  is  along  the  blade  span,  and  which  results  from  temporal 
fluctuations  in  the  loading. 


The  contribution  from  the  axisymmetric  exponential  terms  is  given 
by: 
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which  is  also  proportional  to  <Spo  for  the  reason  cited  above.  and  ?2  are 
defined  by: 


(33) 


The  first  sum  is  over  those  modes  which  are  above  cut-off,  the  last  of  which 
* 

is  denoted  by  k  ;  the  higher  order  modes  in  the  second  sum  are  all  below 

* 

cut-off.  In  general,  k  will  be  a  function  of  the  reduced  frequency,  d3  , 
and  the  operating  parameters  of  the  rotor.  It  will  also  be  different  for 
different  values  of  the  azimuthal  index,  n  (see  pp.  35-36  of  Ref.  2). 
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The  contribution  from  the  asymmetric  exponential  terms  is  given  by: 
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The  significance  of  the  two  sums  and  the  index  k  is  the  same  as  described 

•k 

for  the  axisymmetric  terms,  except  here  of  course  k  =  k(n).  Note  that 

ih  + 

the  algebraic  forms  of  F^,  F?,  G~  and  Q”  remain  the  same  whether  a  mode  is 
above  or  below  cut-off;  all  that  changes  is  the  sign  of  8* ,  . 
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Another  step  that  has  been  taken  in  writing  the  above  results  is 
the  truncation  of  the  azimuthal  and  radial  mode  series.  The  original  limits 
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on  m  and  k  have  been  reduced  from  -  oo  and  co  to  -  NM  and  NK,  respectively. 

The  magnitude  of  these  truncation  limits  directly  affects  the  dimensioning 
of  the  computer  program.  By  estimating  the  magnitude  of  the  radial  eigenvalues,, 
,  and  the  resultant  for  various  extreme  cases,  it  was  felt 

that  roughly  the  same  dimensioning  used  in  the  steady  flow  code  should  be 
sufficient;  viz..  NM  =  5,  corresponding  to  11  azimuthal  modes,  and  NK  =  30 
radial  modes. 

The  full  equation  for  the  normal  velocity  field  is  given  by  the  sum 
of  Eqs.  (28),  (32)  and  (35).  The  next  step  is  the  specification  of  tC  ,  in 
terms  of  the  prescribed  inflow  distortion  pattern,  at  a  discrete  set  of 
collocation  points  on  the  reference  blade  surface  denoted  by: 

<r  =  <r -  n  =  i ,  2  ...  n  N 

X  =  ;  ,  2  , . .  NL  (36) 


At  each  point  two  equations  will  result,  one  each  for  the  real  and  imaginary 
parts  of  VA  .  Thus:,  if  NN  x  NL  is  chosen  equal  to  NI  x  NJ,  a  determinate  set 
of  2  x  NI  x  NJ  simultaneous,  linear  algebraic  equations  will  result.  This 
system  can  then  be  viewed  as  a  single  matrix  equation  to  be  solved  for  the 
column  vector  of  unknows  consisting  of  cl.-  and  h- .  The  rank  of  the  matrix 
may  be  reduced  by  recognizing  that  the  hard-wall  boundary  conditions  require 
that  b  /  her  vanish  at  both  the  hub  and  tip.  This  can  be  used  directly 
in  Eq.  (17)  to  reduce  the  number  of  unknowns  (and  hence  collocation  points)  to 
2  x  NI  x  (NJ-2),  as  was  done  in  Ref.  2.  The  reduced  system  is  then  solved 
using  Gaussian  elimination  with  pivoting.  The  local  loading  then  follows 
directly  from  Eq.  (17).  The  real  and  imaginary  parts  of  the  sectional  lift  and 
moment  coefficients,  CL(<r)  and  C^(cr)r  are  obtained  from  Eqs.  (19)  and  (20) 
of  Ref.  2,  with  a. '  ^  replaced  by  (  +  l  h  ^  )  . 


6.  ACOUSTIC  FARFIELD 

In  addition  to  the  unsteady  aerodynamic  loading,  this  model  is  also  cap¬ 
able  of  predicting  the  acousting  farfield  within  the  duct.  This  field  is  composed 
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of  the  modes  above  cut-off,  which  propagate  undamped  away  from  the  rotor.  We 
will  now  show  how  this  sound  field  may  be  directly  related  to  the  loading 
coefficients,  CL 'L  and  b  •  ^  ,  which  result  from  solving  the  integral  equation. 


The  starting  point  is  Eq.  (51)  of  Ref.  2  which  expresses  the  pressure 
perturbation  in  blade-fixed  coordinates.  This  is  transformed  to  duct-fixed 
coordinates  by  changing  &  to  ,  and  only  those  modes  above  cut-off, 

■Xr 

i.e.  k  <  k  ,  are  retained.  The  result  is 
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where  Q.  is  a  constant  defined  by 

aa’^  =  £n  <fiT  r)r  +  M  ^  ±  +■  )7?r/3nk2  / /32  (37b) 

^  K 

Note  that  the  n  =  omode,  if  present,  has  been  included  in  the  sum,  as  has  the 
k  =  o  radial  mode.  It  should  be  recalled,  however,  that  the  latter  is  trivially 
zero  (  kno  ~  0  )  unless  r\  =  o  .  The  superscripts  u,  d  are  used  to  indicate  up- 
or  downstream  radiation,  and  correspond  to  using  the  upper/lower  signs  on  the 
right-hand  side  of  the  equation.  Eq.  (37)  is  just  a  superposition  of  propagating 
duct  acoustic  modes.  All  of  the  information  on  the  spatial  variations  and 
propagation  characteristics  of  the  waves  is  contained  in  the  radial  eigenfunctions 
R,nk  and  in  the  complex  exponential. 

On  the  other  hand,  all  of  the  information  on  how  the  noise  was 
generated  is  contained  in  the  mode  amplitudes,  ,  which  are  independent  of 

both  space  and  time.  The  latter  are  complex,  and  are  expressible  as  an  integral 
of  the  loading  over  the  rotor  blade  surfaces.  The  form  of  the  loading  expansion 
series  again  allows  the  radial  and  chordwise  integrations  to  be  done  analytically. 
Without  giving  the  intermediate  details,  one  obtains  the  following  relationship 
between  the  acoustic  mode  amplitudes,  ,  and  the  loading  expansion 

coefficients,  a...  and  b.  .  : 

V  7 
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Here  the  X  ^  are  the  same  radial  integrals  encountered  in  the  integral 
rxK 

equation,  and  so  will  have  already  been  evaluated  in  the  determination  of  the 
blade  loading.  The  axial  integrals,  on  the  other  hand,  are  even  simpler  than 
those  in  the  integral  equation,  because  in  the  farfield  all  of  the  dependence 
on  X  can  be  factored  out  (cf.  Eq.  (37)).  They  are  defined  by 
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where  the  V/.(Xa)  are  defined  in  Eq.  (24).  These  integrals  may  be  done  analyti¬ 
cally  following  the  same  procedures  as  before,  and  the  results  are  quoted  at 
the  end  of  Appendix  A. 


Eqs .  (37) -(39)  allow  the  evaluation  of  the  unsteady  pressure  at  any 
point  in  the  farfield.  In  practice,  however,  the  quantities  of  most  interest  are 
the  mean  square  pressure  at  any  point,  and  the  total  power  radiated  away  from 
the  rotor.  The  former  is  defined  as 


u  d 

<^>  5  T 


T 

J  [6U?*a  ]  dt 


(40) 


For  the  case  of  present  interest,  i.e.,  the  response  to  a  steady  inlet  distortion, 
only  harmonics  of  the  blade  passage  frequency  will  occur,  and  hence,  T -  (&SL)  . 
Substituting  from  Eq.  (37),  one  obtains 
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Note  that  the  right-hand  side  of  Eq.  (41)  still  depends  on  O'  ,  X  ,  and  <9-  , 

as  it  should. 


From  Ref.  27,  the  total  power  radiated  upstream  or  downstream  of  the 
rotor,  r  ,  may  be  expressed  as 
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The  last  term  is  seen  to  be  proportional  to  the  integral  of  Eq.  (41)  over  the 
duct  cross-section.  This  is  easily  performed  by  taking  advantage  of  the  orthog¬ 
onality  properties  of  the  duct  modes  in  the  azimuthal  and  radial  coordinates,  and 
the  contribution  to  Eq.  (42)  is  found  to  be 


TTrjrt  V  t 
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The  other  terms  in  Eq.  (42)  involve  the  perturbation  field  of  the  axial 
velocity  component,  V"£  .  Since  2TZ  satisfies  the  same  differential  equation 
as  -0  ,  it  too  can  be  expressed  as  a  superposition  of  duct  modes  of  the  same 

'  u  ct 

form  used  in  Eq.  (37).  The  mode  amplitudes  for  iTz  waves,  say  IT  1  ,  can  be 

related  to  the  pressure  mode  amplitudes  through  the  axial  momentum  equation 
(Ref.  28,  pp.  30-31),  viz., 


where 
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is  the  dimensionless  acoustic  admittance  for  the  (  rt)  k  )  mode. 
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With  Eqs.  (43)  and  (44)  one  can  immediately  write  down  an  expression  for 
analogous  to  Eq.  (37)  for  the  pressure.  Expressions  for  and  V’J'  re¬ 

quired  in  the  first  two  terms  of  Eq.  (42)  are  then  formed.  The  integrations 
over  time  and  the  duct  cross-section  go  through  just  as  before,  since  the  only 

.  u  cL 

additional  factors  are  the  A  ’  ,  which  are  independent  of  space  and  time. 
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Finally  then,  the  full  result  for  r  may  be  written  as 
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(45) 


54 


We  note  here  that  Eqs.  (41)  and  (45)  agree  with  our  previous  results  in 
the  rotor-stator  interaction  problem  (Ref.  28,  Eqs.  (22)  and  (28)),  if  there  one 
sets  the  mode  amplitudes  generated  by  one  of  the  two  rows  equal  to  zero. 

7 .  CODE  DEVELOPMENT 

In  developing  the  associated  computer  program  on  the  AFWAL/ASD  computers, 
a  storage  limitation  was  encountered.  As  a  result,  the  maximum  allowable  values 
for  each  of  NL  and  NN  (the  number  of  chordwise  and  spanwise  collocation  points), 
and  NI  and  NJ  (the  number  of  chordwise  and  spanwise  terms  in  the  loading  expan¬ 
sion)  ,  were  reduced  from  10  down  to  5.  This  lowered  the  total  program  storage 
requirements  down  to  about  76K  words.  Even  at  that,  and  despite  the  fact  that 
CPU  time  on  the  ASD  CYBER  74  was  typically  less  than  a  minute,  turnaround  time 
was  limited  to  about  one  run  a  day.  This  could  prove  to  be  a  problem  in  applying 
the  code  to  high  reduced  frequencies,  where  more  collocation  points  and  loading 
expansion  terms  would  likely  be  required. 

Another  issue  which  had  to  be  addressed  concerned  the  axial  integrals, 
of  which  there  are  several  thousand,  due  to  the  various  combinations  of  i,  1,  n, 
and  k  which  can  occur.  One  of  the  principal  advantages  of  our  approach  is  that 
we  have  been  able  to  derive  exact  analytical  expressions  for  these  integrals 
(Appendix  A)  which  results  in  an  enormous  decrease  in  the  amount  of  time  needed 
for  their  evaluation.  Unfortunately,  comparisons  with  numerical  evaluations  of 
the  same  integrals  showed  that  a  portion  of  them  was  being  calculated  inaccurately 
due  to  round-off  error  in  the  arithmetic.  The  problem  is  that  the  integrals 
naturally  express  themselves  as  a  sum  of  modified  Bessel  functions.  As  n  and  k 
increase  in  magnitude,  the  magnitude  of  the  integrals  should  monotonically 
decrease.  But  the  size  of  the  individual  terms  in  the  sum  grows  exponentially, 
and  when  the  disparity  between  the  two  approaches  14  significant  digits  (CDC 
single  precision)  all  accuracy  is  lost.  Essentially,  we  find  ourselves  trying 
to  take  very  small  differences  of  very  large  numbers.  A  similar  difficulty  was 
encountered  in  the  earlier  steady  flow  code  (see  Ref.  2,  Appendix  C) . 
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Accordingly,,  an  effort  was  made  to  develop  an  asymptotic  evaluation  of 
these  integrals  which  would  be  valid  in  the  limit  where  the  appropriate  para¬ 
meter  becomes  large.  For  those  integrals  with  i  >  Z  ,  an  asymptotic  expansion 
was  found  which  yielded  acceptable  accuracy,  i.e.,  about  1%  or  better.  For  i  =  1  , 

the  leading  edge  singularity  was  first  subtracted  out  and  integrated  analytically, 
with  the  asymptotic  analysis  applied  to  the  remaining  integral.  Unfortunately, 
the  resulting  combined  expression  did  not  provide  the  desired  1%  accuracy. 

Further,  it  does  not  help  to  carry  out  the  asymptotic  expansion  to  more  terms, 
since,  as  is  often  the  case  with  such  series,  it  does  not  converge. 

Because  of  this  lack  of  an  overlap  region  for  i  =  1,  where  both  the 

analytical  and  asymptotic  methods  could  be  expected  to  agree,  and  the  need  to 

expedite  program  development,  it  was  decided  to  use  a  numerical  quadrature  scheme 

for  all  those  integrals  where  round-off  error  would  compromise  the  analytical 

results.  An  efficient  scheme,  which  minimizes  the  number  of  integrand  evaluations 

29 

required  to  achieve  a  specified  accuracy,  was  chosen  for  this  purpose.  This 
increased  the  program  running  time  somewhat,  but  not  as  much  as  had  been  expected. 
In  any  case,  run  times  were  still  well  below  what  one  would  expect  a  finite- 
difference  solution  to  the  3-D  equations  to  take, 


8.  NUMERICAL  RESULTS 


Our  choice  for  an  inlet  distortion  pattern  was  one  which  is  steady  in 
duct  coordinates,  independent  of  X  and  o-,  and  consists  solely  of  a  sinusoidal 
variation  in  9  of  the  axial  perturbation  velocity,  v  .  The  number  of  cycles  in 
the  variation  is  denoted  by  n  ,  and  its  amplitude  by  (v,  /  (j)  =  £  «  1. 

Such  a  variation  can  be  viewed  as  a  single  Fourier  component  of  a  more  complex 
distortion  pattern.  If  such  a  pattern  is  expressed  in  blade-fixed  coordinates, 
the  distortion  component  felt  by  the  blade  normal  to  the  undisturbed  stream 
direction,  v  ,  can  be  shown  to  be: 


=  £  Um 


±r°L 


■b  (<t>T<r)‘ 


-  i  <3 


(cut  —  CO  S  } 


(46) 


where  S  =  (S  -  c/2,) / (e/2, ^  =  x  is  the  non-dimensional  chordwise  coordinate. 

The  reduced  frequency  is  related  to  r\0  by  co  -  n  (f>  yT  .  The  radially 
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dependent  factor  in  brackets  simply  reflects  the  fact  that  the  orientation  of  ‘ 
if  and  Ja  relative  to  f  and  U  varies  with  radius. 

No  exact  solutions  are  available  for  the  response  of  a  full  3-D  blade 

row  to  such  a  pattern.  However,  at  high  hub/tip  ratio,  the  results  should  approach 

those  available  for  2-D  flows.  The  simplest  such  flow,  attained  in,  the  limit  of 

low  solidity  and  low  Mach  number,  is  the  response  of  an  isolated  airfoil  in 

incompressible  flow  to  a  convected  gust  of  the  form  given  by  Eq.  (46)  .  The 

30 

dynamic  response  in  this  case  has  been  expressed  in  closed  form  by  Sears. 

For  the  gust  amplitude  of  Eq.  (46),  a  stripwise  application  of  Sears  analysis 
yields  for  the  lift  coefficient, 

q  —  L'D/dfi.U'C') 

_  _  _  (47) 

=  2  rr  <pT  <r  £  «/  /  *  (<tr<r)z  c~iu>S(c3) 
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where  0(u>)  is  the  so-called  Sears  function,  which  has  been  tabulated  by  Kemp. 

In  order  to  compare  the  present  program's  predictions  against  Eq.  (47), 
it  was  run  with  the  following  inputs:  M  =  0.01,  h  =  0.8,  B  =•  30,  <fiT  =  1.0, 

-qT  =  0.05  n^  =  10  and  £  =  0.01.  The  conditions  correspond  to  a  reduced  frequency 
£  =  0.5.  In  Eq.  (17),  NI  and  NJ  were  each  5,  and  NL  =  5 ,  NN  =  3  collocation 
points  were  used.  The  azimuthal  and  radial  mode  series  were  truncated  at 
NM  =  5  and  NK  =  20,  respectively.  The  two  calculations  are  compared  in  Fig.  24a, 
which  shows  the  magnitude  of  the  complex  lift  coefficient  as  a  function  of  span. 
Very  good  agreement  is  exhibited  right  at  the  hub,  but  the  two  predictions 
diverge  significantly  as  the  tip  region  is  approached. 

The  chordwise  load  distribution  in  such  a  flow  is  just  the  same  flat- 
plate  loading  exhibited  in  steady  flow  -  i.e.,  the  entire  chord  responds  in-phase. 
This  distribution  is  compared  with  that  predicted  by  the  present  code  at  the 
hub  in  Fig.  24b.  The  agreement  is  seen  to  be  excellent.  Unfortunately,  the 
results  for  the  phase  of  showed  poor  agreement:  the  phase  given  by  Eq.  (47) 
is  independent  of  radius  and  equal  to  -33.4°,  while  that  predicted  by  the  code 
varied  monotonically  from  +129.9°  at  the  hub  to  +134.4°  at  the  tip. 
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In  attempting  to  resolve  the  discrepancies  in  the  predictions,  it 
was  realized  that,  again  based  on  strip  theory,  the  radial  variation  of  the 
gust  amplitude  in  Eq.  (46)  would  result  in  a  nearly  linear  spanwise  variation 
in  the  circulation,  P  .  This  would  result  in  a  trailing  vortex  pattern  of 
uniform  strength,  which  is  not  accounted  for  by  strip  theory,  but  whose  influence 
is  present  in  the  3-D  code.  In  an  attempt  to  eliminate  this  effect,  the  same 
case  was  re-run  with  the  amplitude  in  Eq.  (46)  (and  hence,  also  (47))  multiplied 
by  (  0  C7”  )  The  resulting  comparison  of  / CLl  is  shown  in  Fig.  25.  As 
expected,  both  calculations  display  less  spanwise  variation  than  in  Fig.  24; 
however,  now  the  strip  theory  predictions  are  8-15%  higher  than  the  present 
results  over  the  entire  span.  The  phase,  which  in  the  strip  theory  is  unaffected 
by  the  change,  is  now  predicted  by  the  code  to  vary  from  131.3°  to  133.3°. 

Efforts  to  improve  the  agreement  by  increasing  the  hub/tip  ratio  and/or 
lowering  the  reduced  frequency  were  unsuccessful.  It  may  be  that  a  programming 
error  remains  in  the  code,  or  perhaps  an  insufficient  number  of  azimuthal  and 
radial  modes  were  included  in  the  calculation.  At  present,  we  believe  the  latter 
is  more  likely,  but  to  correct  it  would  require  further  increases  in  computer 
storage  requirements.  For  the  reasons  outlined  in  Section  IV.  G.,  this  was  not 
practicable  on  the  present  system.  Accordingly,  the  remaining  effort  was  con¬ 
centrated  on  the  development  of  the  finite  difference  code  described  in  Section  B. 
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SECTION  V 

SUMMARY  AND  CONCLUSIONS 


This  report  presents  the  results  of  a  research  program  on  rotating 
stall  in  axial  flow  compressors.  Three  major  tasks  are  described  which  were: 

(1)  the  development  of  an  Euler  code  for  rotor  rows,  (2)  the  development  of  a 
two-dimensional  stability  theory  for  the  compressible  flow  through  two  blade 
rows,  and  (3)  the  development  of  an  unsteady  lifting  surface  theory  for  the 
flow  through  a  three-dimensional  blade  row. 

The  two-dimensional  version  of  the  Euler  code  developed  works  well  over 
a  range  of  subsonic  inlet  conditions  including  cases  where  shock  waves  form  in 
the  blade  passages.  The  code  works  with  either  a  Kutta  condition  or  conditions 
at  downstream  infinity  specified.  The  integration  code  requires  approximately 

-3 

2.5  x  10  seconds  per  grid  point  per  time  step  on  a  Cyber  750  computing  machine 
and  required  160,000  octal  words  of  central  memory.  For  the  grid  used  which  was 
10  by  40  approximately,  800  to  1200  time  steps  were  required  to  reach  the  steady 
state.  The  cases  with  shock  waves  required  the  longer  times.  The  size  of  the 
time  steps  used  correspond  to  Courant  numbers  of  5  to  10.  The  upper  limit  for 
convergence  was  not  established;  however,  the  results  for  an  inlet  Mach  number 
of  .76  appeared  to  suffer  excessive  rothalpy  loss  near  the  leading  edge  with  the 
larger  time  step  sizes.  The  results  for  an  inlet  Mach  number  of  .5  did  not  show 
such  effects  so  that  computing  times  for  these  cases  can  probably  be  reduced 
substantially. 

Although  the  basic  integration  methods  used  for  the  Euler  equations 
have  been  well  publicized  in  the  external  aerodynamics  literature,  several  modi¬ 
fications  were  found  to  be  crucial  in  order  to  apply  the  method  to  internal 
flows.  The  calculation  of  the  metrics  of  the  transformation  proves  to  be 
extremely  important  and  a  revision  of  the  numerical  viscosity  treatment  enlarges 
and  enhances  the  domain  where  converged  solutions  can  be  obtained.  In  particular, 
it  was  found  that  the  metrics  must  be  discretized  in  the  same  spatial  fashion  as 
the  governing  partial  differential  equation  in  order  to  avoid  introducing  source¬ 
like  terms  which  would  quickly  destroy  the  solution. 
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A  small  disturbance  stability  analysis  has  been  presented  for  the  two- 
dimensional,  compressible  flow  through  a  rotor-stator  combination.  The  theory 
has  been  implemented  into  a  computer  code  which  predicts  the  stability  boundary 
and  propagation  speed  of  disturbances  for  given  steady  state  blade  row  performance 
data.  The  theory  and  numerical  results  reduce  to  the  incompressible  flow  case 
as  inlet  Mach  number  approaches  zero.  The  calculations  show  that  even  for  fixed 
blade  row  performance,  compressibility  has  a  significant  effect  on  the  stability 
boundary.  It  is  generally  found  that  this  effect  is  destabilizing  compared  to 
the  incompressible  flow  results. 

The  unsteady  lifting  surface  theory  was  developed  for  the  purpose  of 
predicting  the  aerodynamic  response  of  an  annular  blade  row  to  inlet  flow 
distortion.  A  linearized  but  fully  three-dimensional  lifting  surface  analysis 
was  used.  The  model  is  also  capable  of  predicting  the  acoustic  far  field  within 
the  duct.  Numerical  calculations  were  performed  to  determine  the  unsteady 
spanwise  and  chordwise  load  distribution  on  a  rotor  row  subjected  to  a  circum¬ 
ferentially  sinusoidal  perturbation  of  the  axial  inlet  velocity.  In  an  attempt 
to  validate  the  model,  its  results  were  compared  in  the  limit  of  low  Mach 
number,  low  solidity,  and  high  hub/tip  ratio  against  a  stripwise  application  of 
Sears  analysis  for  an  isolated  airfoil  subjected  to  a  sinusoidal  gust.  The 
predicted  magnitudes  of  the  unsteady  lift  coefficients  agree  reasonably  well 
near  the  hub,  and  the  chordwise  load  distributions  there  are  in  excellent  agree¬ 
ment.  However,  the  discrepancies  between  the  two  predictions  grow  as  the  tip 
region  is  approached.  Also,  there  are  major  differences  in  the  phase  of  the  lift 
coefficient  predicted  by  the  two  methods.  The  reason  for  these  differences  is 
unresolved  at  present. 
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Appendix  A 

EVALUATION  OF  AXIAL  INTEGRALS 

This  appendix  summarizes  the  analytical  expressions  obtained  for 
the  axial  integrals  which  arise  in  the  three-dimensional  unsteady  lifting- 
surface  theory.  In  what  follows,  (p  =  X  ,  and  denotes  the  modified 

Bessel  function  of  order  n  (not  to  be  confused  with  the  azimuthal  mode 
number).  In  general,  In  will  be  complex;  its  argument  will  vary  according 
to  the  class  of  integrals  being  evaluated,  as  will  be  described  below. 

For  the  wake  term  integrals,  Eq.  (25),  one  obtains: 
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where  all  the  I„  have  the  imaginary  argument,  (_-tu>  ).  It  should  be  noted 
that  with  such  an  argument,  1^  is  purely  real/imaginary  for  n  even/odd. 

The  results  for  the  exponential  integrals  above  cut-off,  Eq.  (26), 

are : 
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where  in  the  above  expressions 

In  =  In  Ca.  ±b)  1 

where  a  and  b  are  defined  in  Eq.  (26c) . 

The  results  for  the  exponential  integrals  below  cut-off,  Eq.  (27), 

are: 
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where  the  argument  of  In  in  the  above  expressions  is  (  C-icL'),  where  c  and 
oL  are  defined  in  Eq.  (27c)  . 

In  the  steady  flow  problem,  for  a  2-  2  it  was  found  that  the  axial 
integrals  evaluated  at  -  X  differed  from  those  atf£  by  at  most  a  change  in 
sign,  depending  on  4-  (see  Appendix  C  of  Ref.  2) .  Not  surprisingly, 
these  symmetry  properties  still  hold  true  in  the  present  case  for  the  expo¬ 
nential  integrals  below  cut-off.  For  those  above  cut-off  and  for  the  wake 
term  integrals,  the  relationships  are  not  so  simple.  In  fact,  they  are  no 
longer  "symmetry"  properties  in  the  strict  meaning  of  the  word,  but  for  want 
of  a  better  terminology,  they  will  be  referred  to  here  as  such. 
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For  the  integrals  in  the  wake  terms,  one  finds  that 
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For  the  axial  integrals  pertaining  to  modes  above  cut-off,  one 
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where  again  I ^  indicates  1^  [~  <■  (  &  ±  b  )  ] 

Finally,  for  the  axial  integrals  pertaining  to  modes  below  cut-off, 
we  have  for  i  2.  2  , 
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The  symmetries  described  in  Eqs.  (A- 21)  through  (A-34)  can  yield  a 
big  savings  in  computer  time  if  the  axial  collocation  points  are  symmetrically 
disposed  about  the  mid-chord,  X  =  0  .  Everything  on  the  right-hand  side  of 
these  equations  will  be  known  from  having  evaluated  the  integrals  at  +x 
The  other  half,  at-X  ,  can  then  be  evaluated  using  these  relations  with  very 
little  effort. 

The  analytical  expressions  given  in  this  appendix,  though  admittedly 

complicated,  are  nevertheless  preferable  to  using  numerical  quadrature.  This 
might  not  be  true  if  one  had  only  to  perform  on  the  order  of  ten,  or  perhaps 

even  a  hundred,  such  integrals.  But  when  literally  thousands  are  required, 
as  is  the  case  here,  the  systematic  way  in  which  these  expressions  separate 
the  influence  of  the  indices  n  ,  k  ,  l  and  X  affords  advantages  in  programming 
that  outweigh  their  complexity.  For  example,  once  the  collocation  points  X.£ 
are  chosen,  the  corresponding  values  of  <t>  and  all  the  sines  can  be  determined. 
Similarly,  note  that  the  arguments  of  the  in  the  integrals  arising  from 
the  exponential  terms  depend  only  on  the  azimuthal  and  radial  mode  indices, 
n  and  k  .  Since  the  argument  is  independent  of  the  order  of  the  Bessel 
function,  one  can  make  use  of  a  very  efficient  algorithm  based  on  a  recursion 
relation  (Ref.  32)  to  calculate  all  orders  of  the  function  at  once  for 
a  given  mode.  A  subroutine  to  do  this  had  already  been  written  and  debugged 
as  part  of  our  work  on  the  steady  flow  problem.  Furthermore,  the  evaluation 
of  the  analytical  expressions  renders  the  application  of  the  symmetry  pro¬ 
perties  in  Eqs.  (A-21)  through  (A-26)  almost  trivial,  which  would  not  be 
true  if  numerical  quadrature  were  employed  exclusively.  For  these  reasons 
the  use  of  the  expressions  given  above  yields  a  big  savings  in  computer  time 
over  quadrature  schemes,  which  require  that  each  step  of  the  evaluation  be 


(A- 32 ) 
(A- 33} 
(A-34) 
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repeated  for  all  combinations  of  rt  ,  k  ,  i  ,  and  J.  .  On  the  other  hand, 
round-off  error  can  ruin  the  accuracy  of  these  expressions  for  some  modes; 
this  point  is  discussed  in  the  main  text. 

Finally,  we  quote  the  results  for  the  axial  integrals  which  appear 
in  the  expressions  for  the  acoustic  farfield,  Eqs .  (38) -(39).  These  are: 

i  =  1 


a,  d 


(C)*)'  =  7T  (<•/-) 


(51*)“'  - 


i  >  2  and  even 
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(A- 35b) 


i  >  2  and  odd 


(ci*)  ' 
(s:*)“’ 


=  o 


(-  0  i 
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L-i  '  n* 


*:r  )  +  XL  ( o.unkd ) 


(A- 35  c) 


where  J".  is  the  Bessel  function  of  the  first  kind,  of  order  i. 
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APPENDIX  B 

GOVERNING  INTEGRAL  EQUATION 


For  easy  reference,  we  reproduce  here  the  complete  integral  equation 
relating  the  velocity  component  normal  to  the  blades,  i rn  ,  to  the  unknown 
loading,  : 
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The  remaining  symbols  are  defined  in  the  main  text.  For  a  derivation  of  this 
equation,  the  reader  is  referred  to  Ref.  2. 
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Figure  1  COMPUTATIONAL  DOMAIN 


Figure  2  COMPUTATIONAL  GRID  IN  PHYSICAL  PLANE  FOR  TW-1  CASCADE  SOLIDITY,  a  =  1.16 
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Figure  3  COMPUTATIONAL  GRID  IN  PHYSICAL  PLANE 
NACA  65  (12)  10  CASCADE 
SOLIDITY,  a  =  1.0 
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NORMALIZED  SURFACE  PRESSURE,  P J  p 
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TRAILING  UPPER  £  INDEX  LOWER  TRAILING 
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PRESSURE  COEFFICIENT,  C„  =  -» — ~  cos 
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Figure  5  SURFACE  PRESSURE  COEFFICIENT  FOR  TW-1  CASCADE 
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NORMALIZED  SURFACE  PRESSURE, 
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Figure  7  SURFACE  PRESSURE  COEFFICIENT  FOR  NACA  65  (12)  10  CASCADE 


PRESSURE  COEFFICIENT,  C 


Figure  8  SURFACE  PRESSURE  DISTRIBUTION  NACA  65  (12)  10  CASCADE 
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PRESSURE  COEFFICIENT,  C 
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PRESSURE  COEFFICIENT,  C 


Figure  11  SURFACE  PRESSURE  DISTRIBUTION  NACA  65  (12)  10  CASCADE 


PRESSURE  COEFFICIENT,  C 


1.0 


Figure  12  SURFACE  PRESSURE  DISTRIBUTION  NACA  65  (12)  10  CASCADE 
WITH  AND  WITHOUT  KUTTA  CONDITION 
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PRESSURE  COEFFICIENT,  C 


Figure  13  SURFACE  PRESSURE  DISTRIBUTION  NACA  65  (12)  10  CASCADE  WITH 
AND  WITHOUT  KUTTA  CONDITION 
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SOLID  LINES  -  PRESENT  CALCULATIONS 
SYMBOLS  •  CALCULATIONS  FROM  REF.  5 


Figure  15  THEORETICAL  STABILITY  CHARACTERISTICS  OF  STATOR  SET  NO.  4, 
SINGLE  BLADE  ROW  THEORY 


PROPAGATION  VELOCITY 


INLET  RELATIVE  SWIRL, 


Figure  16  THEORETICAL  PROPAGATION  VELOCITIES  FOR  STATOR  SET  NO.  4, 
SINGLE  BLADE  ROW  THEORY 
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DAMPING  FACTOR 


Figure  17  THEORETICAL  STABILITY  CHARACTERISTICS  OF  ROTOR  NO.  1, 
SINGLE  BLADE  ROW  THEORY 
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DAMPING  FACTOR 


Figure  19  THEORETICAL  STABILITY  CHARACTERISTICS  OF  ROTOR-STATOR  STAGE 
STATOR  STAGGER  ANGLE,  6SM  =  37.2  deg 
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Figure  20  INFLUENCE  OF  AXIAL  MACH  NO.  ON  STAGE  STABILITY  CHARACTERISTICS 
OF  ROTOR-STATOR  STAGE;  <5  cM  =  28.2°,  n  =  1 


DAMPING  FACTOR,  ^h- 


Figure  21  INFLUENCE  OF  AXIAL  MACH  NO.  ON  STAGE  STABILITY  CHARACTERISTICS 
OF  ROTOR-STATOR  STAGE;  <5  SM  =  37.2°,  n  =  1 
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Figure  25  COMPARISON  OF  PRESENT  THEORY  VERSUS  2-0  STRIP  THEORY  FOR  CASE 
WITH  NEARLY  CONSTANT  CIRCULATION 
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